! $Id: HOME_kernels_cne.F90,v 1.1 2009/10/14 21:47:50 nnz Exp nnz $

#define SV_   sv(2,8)
#define IN_   sv(1,1):sv(1,2)
#define UN_   sv(1,3):sv(1,4)
#define INH_  sv(1,5):sv(1,6)
#define IS_   sv(1,7)
#define IE_   sv(1,8)
#define JN_   sv(2,1):sv(2,2)
#define VN_   sv(2,3):sv(2,4)
#define JNH_  sv(2,5):sv(2,6)
#define JS_   sv(2,7)
#define JE_   sv(2,8)

#define INLO_ sv(1,5)
#define INHI_ sv(1,6)
#define JNLO_ sv(2,5)
#define JNHI_ sv(2,6)

!BOI
module HOME_kernels_cne

!MODULE: HOME_kernels_cne

!DESCRIPTION:
! "HOME_kernels" provides low-level numerical kernels for computing 
! intermediate quantities (fluxes, derivatives, etc...) for use in
! the dynamical core of ocean models with arbitrary vertical coordinates.
!  - Enjoy!

!AUTHOR:
! N.T.Zadeh   (NOAA/GFDL)
! A.J.Adcroft (NOAA/GFDL)

!USES:
! none

!IMPLICIT NONE

contains
!EOI

! ============================================================================

subroutine xTracerFlx_upwind_cne(sv,dy,uhrho,tracer_field,tracer_flux)

!INTERFACE: xTracerFlx_upwind_cne(sv,dy,uhrho,tracer_field,tracer_flux)

!DESCRIPTION: Calculates the upwind tracer flux x direction.
!NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 1 points on each boundary.

!ARGUMENTS:
  integer, intent(in) :: SV_	        	! Shape vector
  real,    intent(in) :: dy(UN_,JN_)    	! y width of face at u-point (m)
  real,    intent(in) :: uhrho(UN_,JN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(UN_,JN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
!LOCAL VARIABLES:
  integer :: i,j
  real :: vp,vm
 
  do j=JS_,JE_; do I=IS_,IE_
    vp = (uhrho(I,j) + ABS(uhrho(I,j))) * 0.5
    vm = uhrho(I,j)-vp
    tracer_flux(I,j) = dy(I,j) * (vp * tracer_field(i,j) + vm * tracer_field(i+1,j))
  enddo; enddo

end subroutine xTracerFlx_upwind_cne

subroutine yTracerFlx_upwind_cne(sv,dx,vhrho,tracer_field,tracer_flux)

!INTERFACE: yTracerFlx_upwind_cne(sv,dx,vhrho,tracer_field,tracer_flux)

!DESCRIPTION: Calculates the upwind tracer flux y direction.
!NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 1 points on each boundary.

!ARGUMENTS:
  integer, intent(in) :: SV_	        	! Shape vector
  real,    intent(in) :: dx(IN_,VN_)    	! x width of face at u-point (m)
  real,    intent(in) :: vhrho(IN_,VN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(IN_,VN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
! Local variables
  integer :: i,j
  real :: vp,vm

  do J=JS_,JE_; do i=IS_,IE_
    vp = (vhrho(i,J) + ABS(vhrho(i,J))) * 0.5
    vm = vhrho(i,J)-vp
    tracer_flux(i,J) = dx(i,J) * (vp * tracer_field(i,j) + vm * tracer_field(i,j+1))
  enddo; enddo

end subroutine yTracerFlx_upwind_cne

! ============================================================================


subroutine xTracerFlx_2ndOrderCentered_cne(sv,dy,uhrho,tracer_field,tracer_flux)

!INTERFACE: xTracerFlx_2ndOrderCentered_cne(sv,dy,uhrho,tracer_field,tracer_flux)

!DESCRIPTION: Calculates the 2nd Order Centered tracer flux in x direction.
!NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 1 points on each boundary.

!ARGUMENTS:
  integer, intent(in) :: SV_		! Shape vector
  real,    intent(in) :: dy(UN_,JN_)    	! y width of face at u-point (m)
  real,    intent(in) :: uhrho(UN_,JN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(UN_,JN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
!LOCAL VARIABLES:
  integer :: i,j

  
  do j=JS_,JE_; do I=IS_,IE_
    tracer_flux(I,j) = dy(I,j) * uhrho(I,j) * (tracer_field(i,j) + tracer_field(i+1,j)) * 0.5
  enddo; enddo

end subroutine xTracerFlx_2ndOrderCentered_cne

subroutine yTracerFlx_2ndOrderCentered_cne(sv,dx,vhrho,tracer_field,tracer_flux)

!INTERFACE: yTracerFlx_2ndOrderCentered_cne(sv,dx,vhrho,tracer_field,tracer_flux)

!DESCRIPTION: Calculates the 2nd Order Centered tracer flux in y direction.
!NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 1 points on each boundary.

!ARGUMENTS:
  integer, intent(in) :: SV_		! Shape vector
  real,    intent(in) :: dx(IN_,VN_)    	! x width of face at u-point (m)
  real,    intent(in) :: vhrho(IN_,VN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(IN_,VN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
!LOCAL VARIABLES:
  integer :: i,j

  
  do J=JS_,JE_; do i=IS_,IE_
    tracer_flux(i,J) = dx(i,J) * vhrho(i,J) * (tracer_field(i,j) + tracer_field(i,j+1)) * 0.5
  enddo; enddo

end subroutine yTracerFlx_2ndOrderCentered_cne


! ============================================================================

! ============================================================================


subroutine xTracerFlx_4thOrderCentered_cne(sv,tmask,dy,uhrho,tracer_field,tracer_flux)

  !INTERFACE: xTracerFlx_4thOrderCentered_cne(sv,tmask,dy,uhrho,tracer_field,tracer_flux)

  !DESCRIPTION: Calculates the 4th Order Centered tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !ARGUMENTS:
  integer, intent(in) :: SV_	        	! Shape vector
  real,    intent(in) :: dy(UN_,JN_)    	! y width of face at u-point (m)
  real,    intent(in) :: uhrho(UN_,JN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tmask(INH_,JNH_)       ! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(UN_,JN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j,im1,ip2
  real, parameter :: a4=7.0/12.0, b4=-1.0/12.0

  do j=JS_,JE_; do I=IS_,IE_
     im1   = int(tmask(i-1,j)*(i-1) + (1.0-tmask(i-1,j))*(i))
     ip2   = int(tmask(i+2,j)*(i+2) + (1.0-tmask(i+2,j))*(i+1))  

     tracer_flux(I,j) = dy(I,j) * uhrho(I,j) * &
          (a4*(tracer_field(i+1,j)+tracer_field(i,j)) + &
           b4*(tracer_field(ip2,j)+tracer_field(im1,j)))

  enddo; enddo

end subroutine xTracerFlx_4thOrderCentered_cne

subroutine yTracerFlx_4thOrderCentered_cne(sv,tmask,dx,vhrho,tracer_field,tracer_flux)

  !INTERFACE: yTracerFlx_4thOrderCentered_cne(sv,tmask,dx,vhrho,tracer_field,tracer_flux)

  !DESCRIPTION: Calculates the 4th Order Centered tracer flux in y direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !ARGUMENTS:
  integer, intent(in) :: SV_	        	! Shape vector
  real,    intent(in) :: dx(IN_,VN_)    	! y width of face at u-point (m)
  real,    intent(in) :: vhrho(IN_,VN_)  	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) :: tmask(INH_,JNH_)       ! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)!Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(IN_,VN_)	! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j,jm1,jp2
  real, parameter :: a4=7.0/12.0, b4=-1.0/12.0

  do J=JS_,JE_; do i=IS_,IE_
     jm1   = int(tmask(i,j-1)*(j-1) + (1.0-tmask(i,j-1))*(j))
     jp2   = int(tmask(i,j+2)*(j+2) + (1.0-tmask(i,j+2))*(j+1))  

     tracer_flux(i,J) = dx(i,J) * vhrho(i,J) * &
          (a4*(tracer_field(i,j+1)+tracer_field(i,j)) + &
           b4*(tracer_field(i,jp2)+tracer_field(i,jm1)))

  enddo; enddo

end subroutine yTracerFlx_4thOrderCentered_cne


! ============================================================================

subroutine xTracerFlx_mdflSupB_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_mdflSupB_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the multi-dimensional flux-limited tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real    :: Rjm, Rj, Rjp, Cr, cfl, massflux

  tracer_flux = 0.0

  do j=JNLO_,JNHI_; do I=IS_,IE_ 

     Rjp = ( tracer_field(i+2,j) - tracer_field(i+1,j) )       &
          *         tmask(i+2,j) *        tmask(i+1,j)
     Rj =  ( tracer_field(i+1,j) - tracer_field( i ,j) )       &
          *         tmask(i+1,j) *        tmask( i ,j)
     Rjm = ( tracer_field( i ,j) - tracer_field(i-1,j) )       &
          *         tmask( i ,j) *        tmask(i-1,j)

     massflux = dy(I,j) * uhrho(I,j)

     cfl = abs( u(I,j) * dt /dx(I,j)) 

     if ( Rj .NE. 0.0) then
        if ( massflux .GT. 0.0) then
           Cr = Rjm / Rj
        else
           Cr = Rjp / Rj
        endif
     else
        if ( massflux .GT. 0.0) then
           Cr = Rjm * 1.0e20
        else
           Cr = Rjp * 1.0e20
        endif
     endif

     Cr = max(0.0, max(min(1.0, 2.0 * Cr), min(2.0, Cr)))

     tracer_flux(I,j) = 0.5 * ( massflux * (tracer_field(i+1,j) + tracer_field(i,j))  &
          - abs(massflux) * ( 1.0 + (cfl - 1 ) * Cr ) * Rj )    &
          * tmask(i,j) * tmask(i+1,j)


  enddo; enddo

end subroutine xTracerFlx_mdflSupB_cne

subroutine yTracerFlx_mdflSupB_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_mdflSupB_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the multi-dimensional flux-limited tracer flux in y direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out):: tracer_flux(IN_,VN_)   ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real    :: Rjm, Rj, Rjp, Cr, cfl, massflux

  do J=JS_,JE_; do i=INLO_,INHI_

     Rjp = ( tracer_field(i,j+2) - tracer_field(i,j+1) )       &
          * tmask(i,j+2) * tmask(i,j+1)
     Rj =  ( tracer_field(i,j+1) - tracer_field( i ,j) )       &
          * tmask(i,j+1) * tmask( i ,j)
     Rjm = ( tracer_field( i ,j) - tracer_field(i,j-1) )       &
          * tmask( i ,j) * tmask(i,j-1)

     massflux = dx(i,J) * vhrho(i,J)

     cfl = abs( v(i,J) * dt /dy(i,J)) 

     if ( Rj .NE. 0.0) then
        if ( massflux .GT. 0.0) then
           Cr = Rjm / Rj
        else
           Cr = Rjp / Rj
        endif
     else
        if ( massflux .GT. 0.0) then
           Cr = Rjm * 1.0e20
        else
           Cr = Rjp * 1.0e20
        endif
     endif

     Cr = max(0.0, max(min(1.0, 2.0 * Cr), min(2.0, Cr)))

     tracer_flux(i,J) = 0.5 * ( massflux * (tracer_field(i,j+1) + tracer_field(i,j))  &
                                - abs(massflux) * ( 1.0 + (cfl - 1 ) * Cr ) * Rj )    &
                            * tmask(i,j) * tmask(i,j+1)


  enddo; enddo

end subroutine yTracerFlx_mdflSupB_cne

! ============================================================================

subroutine xTracerFlx_plm_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_plm_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise linear tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.


  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: a,b	!  slopes
  real    :: slopeTimesWidth, cfl, massflux

  tracer_flux = 0.0

  do j=JS_,JE_; do I=IS_,IE_ 

     massflux = dy(I,j) * uhrho(I,j)

     cfl = abs( u(I,j) * dt /dx(I,j)) 

     if (massflux .GT. 0.0) then
       a = (tracer_field(i+1,j)-tracer_field(i-1,j))*tmask(i+1,j)*tmask(i-1,j)
       b = (tracer_field(i  ,j)-tracer_field(i-1,j))*tmask(i  ,j)*tmask(i-1,j)
     else
       a = (tracer_field(i+2,j)-tracer_field(i  ,j))*tmask(i+2,j)*tmask(i  ,j)
       b = (tracer_field(i+2,j)-tracer_field(i+1,j))*tmask(i+2,j)*tmask(i+1,j)
     endif

     if ( (a * b) .GT. 0.0 ) then
        slopeTimesWidth = sign ( min (abs(a),abs(b)), a )
     else
        slopeTimesWidth = 0.0
     end if

     if(massflux .gt. 0.0) then        
        tracer_flux(I,j) = massflux * (tracer_field(i  ,j) + slopeTimesWidth *(1.0-cfl)/2.0) 
     else
        tracer_flux(I,j) = massflux * (tracer_field(i+1,j) - slopeTimesWidth *(1.0-cfl)/2.0)    
     endif


  enddo; enddo

end subroutine xTracerFlx_plm_cne

subroutine yTracerFlx_plm_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_plm_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise linear tracer flux in y direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference: Durran's book NMWEGFD, p.276

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: a,b	!  slopes
  real    :: slopeTimesWidth, cfl, massflux

  tracer_flux = 0.0

  do J=JS_,JE_; do i=IS_,IE_ 

     massflux = dx(i,J) * vhrho(i,J)
     cfl = abs( v(i,J) * dt /dy(i,J)) 

     if (massflux .GT. 0.0) then
        a = (tracer_field(i,j+1)-tracer_field( i ,j-1))*tmask(i,j+1)*tmask( i ,j-1)
        b = (tracer_field(i  ,j)-tracer_field(i,j-1))*tmask(i  ,j)*tmask(i,j-1)
     else
        a = (tracer_field(i,j+2)-tracer_field( i ,j))*tmask(i,j+2)*tmask( i ,j)
        b = (tracer_field(i,j+2)-tracer_field(i,j+1))*tmask(i,j+2)*tmask(i,j+1)
     endif

     if ( (a * b) .GT. 0.0 ) then
        slopeTimesWidth = sign ( min (abs(a),abs(b)), a )
     else
        slopeTimesWidth = 0.0
     end if

     if(massflux .gt. 0.0) then        
        tracer_flux(i,J) = massflux * (tracer_field(i  ,j) + slopeTimesWidth *(1.0-cfl)/2.0) 
     else
        tracer_flux(i,J) = massflux * (tracer_field(i,j+1) - slopeTimesWidth *(1.0-cfl)/2.0)    
     endif


  enddo; enddo

end subroutine yTracerFlx_plm_cne

subroutine xTracerFlx_laxwen_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_laxwen_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference: Durran's book NMWEGFD, p.276

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: a,b	!  slopes
  real    :: slopeTimesWidth, cfl, massflux

  tracer_flux = 0.0

  do j=JS_,JE_; do I=IS_,IE_ 

     massflux = dy(I,j) * uhrho(I,j)
     cfl = abs( u(I,j) * dt /dx(I,j)) 

     !Durran's algorithm for calculating slopes.
     a = (tracer_field(i+1,j)-tracer_field( i ,j))*tmask(i+1,j)*tmask( i ,j)
     if(massflux .ge. 0.0) then
        b = (tracer_field(i  ,j)-tracer_field(i-1,j))*tmask(i  ,j)*tmask(i-1,j)
     else
        b = (tracer_field(i+2,j)-tracer_field(i+1,j))*tmask(i+2,j)*tmask(i+1,j)
     endif

     if ( (a * b) .GT. 0.0 ) then
        slopeTimesWidth = sign ( min (abs(a),abs(b)), a )
     else
        slopeTimesWidth = 0.0
     end if

     if(massflux .gt. 0.0) then        
        tracer_flux(I,j) = massflux * (tracer_field(i  ,j) + slopeTimesWidth *(1.0-cfl)/2.0) 
     else
        tracer_flux(I,j) = massflux * (tracer_field(i+1,j) - slopeTimesWidth *(1.0-cfl)/2.0)    
     endif


  enddo; enddo

end subroutine xTracerFlx_laxwen_cne

subroutine yTracerFlx_laxwen_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_laxwen_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the  tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference: Durran's book NMWEGFD, p.276

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: a,b	!  slopes
  real    :: slopeTimesWidth, cfl, massflux

  tracer_flux = 0.0

  do J=JS_,JE_; do i=IS_,IE_ 

     massflux = dx(i,J) * vhrho(i,J)
     cfl = abs( v(i,J) * dt /dy(i,J)) 

     !Durran's algorithm for calculating slopes.
     a = (tracer_field(i,j+1)-tracer_field( i ,j))*tmask(i,j+1)*tmask( i ,j)
     if(massflux .ge. 0.0) then
        b = (tracer_field(i  ,j)-tracer_field(i,j-1))*tmask(i  ,j)*tmask(i,j-1)
     else
        b = (tracer_field(i,j+2)-tracer_field(i,j+1))*tmask(i,j+2)*tmask(i,j+1)
     endif

     if ( (a * b) .GT. 0.0 ) then
        slopeTimesWidth = sign ( min (abs(a),abs(b)), a )
     else
        slopeTimesWidth = 0.0
     end if

     if(massflux .gt. 0.0) then        
        tracer_flux(i,J) = massflux * (tracer_field(i  ,j) + slopeTimesWidth *(1.0-cfl)/2.0) 
     else
        tracer_flux(i,J) = massflux * (tracer_field(i,j+1) - slopeTimesWidth *(1.0-cfl)/2.0)    
     endif


  enddo; enddo

end subroutine yTracerFlx_laxwen_cne


! ============================================================================
subroutine xTracerFlx_plmLaurent_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_plmLaurent_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise linear tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 1 points on each boundary.
  
  !Reference: Laurent's plm.F90

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: sigma_l, sigma_c, sigma_r	! left, center and right 
  real    :: slopeTimesWidth, cfl, massflux

  tracer_flux = 0.0

  do j=JS_,JE_; do I=IS_,IE_

     massflux = dy(I,j) * uhrho(I,j)

     cfl = abs( u(I,j) * dt /dx(I,j)) 

     if (massflux .GT. 0.0) then
       sigma_r = (tracer_field(i+1,j)-tracer_field(i  ,j))    *tmask(i+1,j)*tmask(i  ,j)
       sigma_l = (tracer_field(i  ,j)-tracer_field(i-1,j))    *tmask(i  ,j)*tmask(i-1,j)
     else
       sigma_r = (tracer_field(i+2,j)-tracer_field(i+1,j))    *tmask(i+2,j)*tmask(i+1,j)
       sigma_l = (tracer_field(i+1,j)-tracer_field(i  ,j))    *tmask(i+1,j)*tmask(i  ,j)
     endif
     sigma_c = 0.5*(sigma_r+sigma_l)

     if ( (sigma_r * sigma_l) .GT. 0.0 ) then
       slopeTimesWidth = sign ( min( 2.*min(abs(sigma_l),abs(sigma_r)), abs(sigma_c) ), sigma_c )
     else
       slopeTimesWidth = 0.0
     end if

     if (massflux .gt. 0.0) then        
       tracer_flux(I,j) = massflux * (tracer_field(i  ,j) + slopeTimesWidth *(1.0-cfl)/2.0) 
     else
       tracer_flux(I,j) = massflux * (tracer_field(i+1,j) - slopeTimesWidth *(1.0-cfl)/2.0)    
     endif

  enddo; enddo

end subroutine xTracerFlx_plmLaurent_cne

! ============================================================================
subroutine xTracerFlx_plmGOLD_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_plmGOLD_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the modified piecewise linear tracer flux in x direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference: Halberg GOLD code

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of face at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: maxslope	!  slopes
  real    :: cfl, massflux
  real, dimension(IS_-2:IE_+2) :: slopeTimesWidth

  tracer_flux = 0.0

  do j=JS_,JE_

     slopeTimesWidth = 0.0 

     !GOLD's algorithm for calculating slopes.
     do i=IS_,IE_+1         
        if (ABS(tracer_field(i+1,j)-tracer_field(i,j)) < &
             ABS(tracer_field(i,j)-tracer_field(i-1,j))) then
           maxslope = 4.0*(tracer_field(i+1,j)-tracer_field(i,j))
        else
           maxslope = 4.0*(tracer_field(i,j)-tracer_field(i-1,j))
        endif

        if ((tracer_field(i+1,j)-tracer_field(i,j)) * (tracer_field(i,j)-tracer_field(i-1,j)) < 0.0) then
           slopeTimesWidth(i) = 0.0
        elseif (ABS(tracer_field(i+1,j)-tracer_field(i-1,j))<ABS(maxslope)) then
           slopeTimesWidth(i) = tmask(i,j)*tmask(i-1,j) * &
                0.5*(tracer_field(i+1,j)-tracer_field(i-1,j))
        else
           slopeTimesWidth(i) = tmask(i,j)*tmask(i-1,j) * 0.5*maxslope
        endif

     enddo!i

     do I=IS_,IE_

        massflux = dy(I,j) * uhrho(I,j)

        cfl =  u(I,j) * dt /dx(I,j) 
        if(massflux .gt. 0.0) then        
           tracer_flux(I,j) = massflux * (tracer_field(i  ,j) + slopeTimesWidth(i)  *0.5*(1.0-cfl)) 
        else
           tracer_flux(I,j) = massflux * (tracer_field(i+1,j) - slopeTimesWidth(i+1)*0.5*(1.0+cfl))    
        endif

     enddo!i

  enddo !j

end subroutine xTracerFlx_plmGOLD_cne

subroutine yTracerFlx_plmGOLD_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_plmGOLD_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the modified piecewise linear tracer flux in y direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference: Halberg GOLD code

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: maxslope	!  slopes
  real    :: cfl, massflux
  real, dimension(JS_-2:JE_+2) :: slopeTimesWidth

  tracer_flux = 0.0

  do i=IS_,IE_ !Niki: Had to invert the order of i and j loops        

     slopeTimesWidth = 0.0 !Niki: Nnot sure about this! What should tracer_flux(IE_+1,j) be?

     !GOLD's algorithm for calculating slopes.
     do j=JS_,JE_+1
          
        if (ABS(tracer_field(i,j+1)-tracer_field(i,j)) < &
             ABS(tracer_field(i,j)-tracer_field(i,j-1))) then
           maxslope = 4.0*(tracer_field(i,j+1)-tracer_field(i,j))
        else
           maxslope = 4.0*(tracer_field(i,j)-tracer_field(i,j-1))
        endif

        if ((tracer_field(i,j+1)-tracer_field(i,j)) * (tracer_field(i,j)-tracer_field(i,j-1)) < 0.0) then
           slopeTimesWidth(j) = 0.0
        elseif (ABS(tracer_field(i,j+1)-tracer_field(i,j-1))<ABS(maxslope)) then
           slopeTimesWidth(j) = tmask(i,j)*tmask(i,j-1) * &
                0.5*(tracer_field(i,j+1)-tracer_field(i,j-1))
        else
           slopeTimesWidth(j) = tmask(i,j)*tmask(i,j-1) * 0.5*maxslope
        endif

     enddo!i

     do J=JS_,JE_
        
        massflux = dx(i,J) * vhrho(i,J)

        cfl = v(i,J) * dt /dy(i,J)
        if(massflux .gt. 0.0) then        
           tracer_flux(i,J) = massflux * (tracer_field(i  ,j) + slopeTimesWidth(j)  *0.5*(1.0-cfl)) 
        else
           tracer_flux(i,J) = massflux * (tracer_field(i,j+1) - slopeTimesWidth(j+1)*0.5*(1.0+cfl))    
        endif

     enddo!i

  enddo !j

end subroutine yTracerFlx_plmGOLD_cne
! ============================================================================

subroutine xTracerFlx_ppm_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_ppm_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise parabolic tracer flux in x direction.
  !             The choice of the limiter is as follows:
  !             ppm_hlimiter=0 : No limiter 
  !             ppm_hlimiter=1 : ppm_limit_cw84
  !             ppm_hlimiter=2 : ppm_limit_ifc
  !             ppm_hlimiter=3 : ppm_limit_sh

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  integer, intent(in) :: ppm_hlimiter           ! Choice of the limiter function
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: cfl, massflux
  real    :: Sim2,Sim1,Si,Sip1,Sip2
  real    :: da2,da4,da3m,da3p, dMx, dMn
  real    :: mskm2,mskm1,msk0,mskp1,mskp2
  real,parameter            :: oneSixth=1./6., r24=1./24., r12=1./12.
  real,parameter            :: fourThirds=4./3., twoThirds=2./3.
  real,dimension(INH_,JNH_) :: da, aL, aR, a6, d1m, d1p, d1mm, d1pp

  tracer_flux = 0.0
  da   = 0.0
  d1m  = 0.0
  d1p  = 0.0
  d1mm = 0.0
  d1pp = 0.0
  aL   = 0.0
  aR   = 0.0
  a6   = 0.0

  do j=JNLO_,JNHI_
     do I=IS_-1,IE_+2 

        ! This block calculates the slope in each cell (aka PLM but
        ! with a fourth-order estimate)
        Sip2 = tracer_field(i+2,j)
        Sip1 = tracer_field(i+1,j)
        Si   = tracer_field(i,j)
        Sim1 = tracer_field(i-1,j)
        Sim2 = tracer_field(i-2,j)
        ! Simple unmasked 4th order
        da4 = r12 * ( 8. * ( Sip1 - Sim1 ) + ( Sim2 - Sip2 ) )
        ! Simple 2nd order
        da2 = 0.5 * ( Sip1 - Sim1 )
        ! Simple 3rd order, biased to left
        da3m = r12 * ( 2. * Sim2 - 12. * Sim1 + 6. * Si + 4. *Sip1 )
        ! Simple 3rd order, biased to right
        da3p = r12 * ( - 4. * Sim1 - 6. * Si + 12. *Sip1  - 2. * Sip2 )
        ! Estimate of slope: da(i,j) is the twice the "mismatch" as defined
        ! by Lin which itself is only half of the slope. The factor of two
        ! in the mismatch is confusing so we use simple estimate of slope.
        ! Note: The temperature masks used here are proxies for the advective
        ! flow masks that would be used if they existed. The only functional
        ! difference is that this scheme will not work with thin walls. 
        mskm2 = tmask(i-2,j)
        mskm1 = tmask(i-1,j)
        msk0  = tmask(i,j)
        mskp1 = tmask(i+1,j)
        mskp2 = tmask(i+2,j)
        da(i,j) = ( (   mskm2 * ( 1. - 0.5 * mskp2 ) * da3m       &
             + mskp2 * ( 1. - 0.5 * mskm2 ) * da3p )     &
             + ( 1. - mskm2 ) * ( 1. - mskp2 ) * da2 )
!!! da(i,j) = da4 ! *** UNMASKED ***
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dMx = max( Sip1, Sim1,  Si ) - Si
        dMn = Si - min( Sip1, Sim1, Si )
        da(i,j) = sign(1.,da(i,j)) * min( abs(da(i,j)) , 2. * min( dMx , dMn ) )  &
             * mskm1 * mskp1 * msk0
        ! Calculate the second difference for use in the Huynh limiter later
        d1m(i,j)  = ( Si - Sim1 ) * mskm1
        d1p(i,j)  = ( Sip1 - Si ) * mskp1
        d1mm(i,j) = ( Sim1 - Sim2 ) * mskm1 * mskm2
        d1pp(i,j) = ( Sip2 - Sip1 ) * mskp1 * mskp2
     enddo !i

     do I=IS_,IE_+1
        ! This block estimates the left and right values on cell boundaries
        ! in i-direction
        Si   = tracer_field(i,j)
        Sim1 = Si - d1m(i,j)
        Sip1 = Si + d1p(i,j)
        !       mskm1 = tmask(i-1,j)
        !       mskp1 = tmask(i+1,j)
        ! Left edge: Eq. B2 in Lin 1994, MWR (132)
        aL(i,j) = 0.5 * ( Sim1 + Si ) + oneSixth * ( da(i-1,j) - da(i,j) )        
        ! ... and masks
        !       aL(i,j) = mskm1 * aL(i,j) + ( 1. - mskm1 ) * Si
        ! Right edge: Eq. B2 in Lin 1994, MWR (132)
        aR(i,j) = 0.5 * ( Si + Sip1 ) + oneSixth * ( da(i,j) - da(i+1,j) )        
        ! ... and masks
        !       aR(i,j) = mskp1 * aR(i,j) + ( 1. - mskp1 ) * Si
     enddo !i
  enddo !j

  ! Limit the edge values, aL and aR, for monotonicity
      if (ppm_hlimiter.eq.1) then
        call ppm_limit_cw84(sv, tracer_field, aL, aR)
      elseif (ppm_hlimiter.eq.2) then
       call ppm_limit_ifc(sv, tracer_field, da, aL, aR)
      elseif (ppm_hlimiter.eq.3) then
        call ppm_limit_sh(sv, tracer_field, &
                         d1m, d1p, d1mm, d1pp, aL, aR)
      endif

  do j=JNLO_,JNHI_
     ! NOTE TO SELF(AJA): SHOULD BE ABLE TO ELIMINATE A6
     do I=IS_,IE_+1
        ! Curvature in i-direction                            
        a6(i,j) = 6. * tracer_field(i,j) - 3. * ( aR(i,j) + aL(i,j) )        
     enddo !i
  enddo !j

  do j=JNLO_,JNHI_
     do I=IS_,IE_

        massflux = dy(I,j) * uhrho(I,j)

        cfl = abs( u(I,j) * dt /dx(I,j)) 

        if (massflux.ge.0.0) then                              
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                * ( aR(i,j) + 0.5 * cfl * ( ( aL(i,j) - aR(i,j) )          &
                + ( 1. - twoThirds * cfl ) * a6(i,j) ) )   
        else                                                  
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                * ( aL(i+1,j) + 0.5 * cfl * ( ( aR(i+1,j) - aL(i+1,j) )    &
                + ( 1. - twoThirds * cfl ) * a6(i+1,j) ) ) 
        endif

     enddo !i                                                
  enddo !j



end subroutine xTracerFlx_ppm_cne

subroutine yTracerFlx_ppm_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_ppm_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise parabolic tracer flux in y direction.

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  integer, intent(in) :: ppm_hlimiter           ! Choice of the limiter function
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: cfl, massflux
  real    :: Sim2,Sim1,Si,Sip1,Sip2
  real    :: da2,da4,da3m,da3p, dMx, dMn
  real    :: mskm2,mskm1,msk0,mskp1,mskp2
  real,parameter                           :: oneSixth=1./6., r24=1./24., r12=1./12.
  real,parameter                           :: fourThirds=4./3., twoThirds=2./3.
  real,dimension(INH_,JNH_)  :: da, aL, aR, a6, d1m, d1p, d1mm, d1pp

  tracer_flux = 0.0
  da   = 0.0
  d1m  = 0.0
  d1p  = 0.0
  d1mm = 0.0
  d1pp = 0.0
  aL   = 0.0
  aR   = 0.0
  a6   = 0.0

  do J=JS_-1,JE_+2
     do i=INLO_,INHI_ 

        ! This block calculates the slope in each cell (aka PLM but
        ! with a fourth-order estimate)
        Sip2 = tracer_field(i,j+2)
        Sip1 = tracer_field(i,j+1)
        Si   = tracer_field(i,j)
        Sim1 = tracer_field(i,j-1)
        Sim2 = tracer_field(i,j-2)
        ! Simple unmasked 4th order
        da4 = r12 * ( 8. * ( Sip1 - Sim1 ) + ( Sim2 - Sip2 ) )
        ! Simple 2nd order
        da2 = 0.5 * ( Sip1 - Sim1 )
        ! Simple 3rd order, biased to left
        da3m = r12 * ( 2. * Sim2 - 12. * Sim1 + 6. * Si + 4. *Sip1 )
        ! Simple 3rd order, biased to right
        da3p = r12 * ( - 4. * Sim1 - 6. * Si + 12. *Sip1  - 2. * Sip2 )
        ! Estimate of slope: da(i,j) is the twice the "mismatch" as defined
        ! by Lin which itself is only half of the slope. The factor of two
        ! in the mismatch is confusing so we use simple estimate of slope.
        ! Note: The temperature masks used here are proxies for the advective
        ! flow masks that would be used if they existed. The only functional
        ! difference is that this scheme will not work with thin walls. 
        mskm2 = tmask(i,j-2)
        mskm1 = tmask(i,j-1)
        msk0  = tmask(i,j)
        mskp1 = tmask(i,j+1)
        mskp2 = tmask(i,j+2)
        da(i,j) = ( (   mskm2 * ( 1. - 0.5 * mskp2 ) * da3m       &
             + mskp2 * ( 1. - 0.5 * mskm2 ) * da3p )     &
             + ( 1. - mskm2 ) * ( 1. - mskp2 ) * da2 )
!!! da(i,j) = da4 ! *** UNMASKED ***
        ! Monotonic constraint, see Eq. B2 in Lin 1994, MWR (132)
        dMx = max( Sip1, Sim1,  Si ) - Si
        dMn = Si - min( Sip1, Sim1, Si )
        da(i,j) = sign(1.,da(i,j)) * min( abs(da(i,j)) , 2. * min( dMx , dMn ) )  &
             * mskm1 * mskp1 * msk0
        ! Calculate the second difference for use in the Huynh limiter later
        d1m(i,j)  = ( Si - Sim1 ) * mskm1
        d1p(i,j)  = ( Sip1 - Si ) * mskp1
        d1mm(i,j) = ( Sim1 - Sim2 ) * mskm1 * mskm2
        d1pp(i,j) = ( Sip2 - Sip1 ) * mskp1 * mskp2
     enddo !i
  enddo !j

  do J=JS_,JE_+1
     do i=INLO_,INHI_
        ! This block estimates the left and right values on cell boundaries
        ! in i-direction
        Si   = tracer_field(i,j)
        Sim1 = Si - d1m(i,j)
        Sip1 = Si + d1p(i,j)
        !       mskm1 = tmask(i,j-1)
        !       mskp1 = tmask(i,j+1)
        ! Left edge: Eq. B2 in Lin 1994, MWR (132)
        aL(i,j) = 0.5 * ( Sim1 + Si ) + oneSixth * ( da(i,j-1) - da(i,j) )        
        ! ... and masks
        !       aL(i,j) = mskm1 * aL(i,j) + ( 1. - mskm1 ) * Si
        ! Right edge: Eq. B2 in Lin 1994, MWR (132)
        aR(i,j) = 0.5 * ( Si + Sip1 ) + oneSixth * ( da(i,j) - da(i,j+1) )        
        ! ... and masks
        !       aR(i,j) = mskp1 * aR(i,j) + ( 1. - mskp1 ) * Si
     enddo !i
  enddo !j

  ! Limit the edge values, aL and aR, for monotonicity
      if (ppm_hlimiter.eq.1) then
        call ppm_limit_cw84(sv, tracer_field, aL, aR)
      elseif (ppm_hlimiter.eq.2) then
       call ppm_limit_ifc(sv, tracer_field, da, aL, aR)
      elseif (ppm_hlimiter.eq.3) then
        call ppm_limit_sh(sv, tracer_field, &
                         d1m, d1p, d1mm, d1pp, aL, aR)
      endif

      ! NOTE TO SELF(AJA): SHOULD BE ABLE TO ELIMINATE A6
  do J=JS_,JE_+1
     do i=INLO_,INHI_
        ! Curvature in i-direction                            
        a6(i,j) = 6. * tracer_field(i,j) - 3. * ( aR(i,j) + aL(i,j) )        
     enddo !i
  enddo !j

  do J=JS_,JE_
     do i=INLO_,INHI_                                         

        massflux = dx(i,J) * vhrho(i,J)

        cfl = abs( v(i,J) * dt /dy(i,J)) 

        if (massflux.ge.0.0) then                              
           tracer_flux(i,J) = massflux * tmask(i,j) * tmask(i,j+1) &
                * ( aR(i,j) + 0.5 * cfl * ( ( aL(i,j) - aR(i,j) )          &
                + ( 1. - twoThirds * cfl ) * a6(i,j) ) )   
        else                                                  
           tracer_flux(i,J) = massflux * tmask(i,j) * tmask(i,j+1) &
                * ( aL(i,j+1) + 0.5 * cfl * ( ( aR(i,j+1) - aL(i,j+1) )    &
                + ( 1. - twoThirds * cfl ) * a6(i,j+1) ) ) 
        endif

     enddo !i                                                
  enddo !j


end subroutine yTracerFlx_ppm_cne

!#######################################################################
! <SUBROUTINE NAME="ppm_limit_cw84">
!
! <DESCRIPTION>
!
! Kernel to limit the edge values for PPM following Colella and Woodward, 1984
! 
! Coded by Alistair.Adcroft@noaa.gov
! Apr 2006
!
! </DESCRIPTION>
!
subroutine ppm_limit_cw84(sv, tracer, aL, aR)
implicit none
! Arguments
integer, intent(in) :: SV_            	! Shape vector
real, dimension(INH_,JNH_), intent(in)    :: tracer
real, dimension(INH_,JNH_), intent(inout) :: aL, aR
! Local
real    :: da2,da3m,da3p,da4,Si
integer :: i,j

    do j=JS_-1,JE_+1
      do i=IS_-1,IE_+1
        ! This block monotonizes the parabola by adjusting the left and right values
        ! Limiter from Colella and Woodward, 1984, Eq. 1.10
        Si   = tracer(i,j)
        if ( ( aR(i,j) - Si ) * ( Si - aL(i,j) ) .le. 0. ) then
           aL(i,j) = Si
           aR(i,j) = Si
        endif
        da2 = aR(i,j) - aL(i,j)            ! Difference of edge values
        da4 = 0.5 * ( aR(i,j) + aL(i,j) )  ! Mean of edge values
        da3m = 6. * da2 * ( Si - da4 )
        da3p = da2 * da2
        if ( da3m .gt.  da3p ) aL(i,j) = 3. * Si - 2. * aR(i,j)
        if ( da3m .lt. -da3p ) aR(i,j) = 3. * Si - 2. * aL(i,j)
      enddo !i
    enddo !j

end subroutine ppm_limit_cw84
! </SUBROUTINE> NAME="ppm_limit_cw84"

!#######################################################################
! <SUBROUTINE NAME="ppm_limit_ifc">
!
! <DESCRIPTION>
!
! Kernel to limit the edge values for PPM using the Improved Full Constraint
! (IFC) of Lin, 2004.
! 
! Coded by Alistair.Adcroft@noaa.gov
! Apr 2006
!
! </DESCRIPTION>
!
subroutine ppm_limit_ifc(sv, tracer, da, aL, aR)
implicit none
! Arguments
integer, intent(in) :: SV_            	! Shape vector
real, dimension(INH_,JNH_), intent(in)    :: tracer
real, dimension(INH_,JNH_), intent(inout) :: da, aL, aR
! Local
real    :: Si,ada,sda
integer :: i,j

    do j=JS_-1,JE_+1
      do i=IS_-1,IE_+1
        ! This block monotonizes the parabola by adjusting the left and right values
        ! Improved full constraint (IFC) limiter from Lin, 2004
        Si  = tracer(i,j)
        ada = abs(da(i,j))
        sda = sign(1., da(i,j))
        aL(i,j) = Si - sda * min( ada, abs(aL(i,j)-Si) )
        aR(i,j) = Si + sda * min( ada, abs(aR(i,j)-Si) )
      enddo !i
    enddo !j

end subroutine ppm_limit_ifc
! </SUBROUTINE> NAME="ppm_limit_ifc"

!#######################################################################
! <SUBROUTINE NAME="ppm_limit_sh">
!
! <DESCRIPTION>
!
! Kernel to limit the edge values for PPM following the monotonicity-
! preserving approach of Suresh and Huynh, 1997.
! 
! Coded by Alistair.Adcroft@noaa.gov
! Apr 2006
!
! </DESCRIPTION>
!
! <NOTE>
! About efficiency: ppm_limit_sh() is not as efficient as it would be if
! we wrote a s/r for each direction. The pre-calculation of d1m, d1p, d1mm and
! d1pp duplicates operations and would be much faster if d1m was replaced
! by d1p(i+1). However, in order to re-use this limiter for the all directions
! (to simplify debugging) I have opted for the less efficient form for now. - AJA
! </NOTE>
subroutine ppm_limit_sh(sv, tracer, d1m, d1p, d1mm, d1pp, aL, aR)
implicit none
! Arguments
integer, intent(in) :: SV_            	! Shape vector
real, dimension(INH_,JNH_), intent(in)    :: tracer, d1m, d1p, d1mm, d1pp
real, dimension(INH_,JNH_), intent(inout) :: aL, aR
! Local
real    :: Si,Sim1,Sip1
real    :: x,y,z,w,dM4m,dM4p,qAV,qMP,qUL,qLC,qMD,qMin,qMax
real, parameter :: fourThirds = 4./3.
integer :: i,j

    do j=JS_-1,JE_+1
      do i=IS_-1,IE_+1
        ! This block monotonizes the parabola by adjusting the left and right values
        ! Limiter from Suresh and Huynh, 1997
        Si   = tracer(i,j)
!       d1m = ( Si - Sim1 ) * tmask(i-1,j)
!       d1p = ( Sip1 - Si ) * tmask(i+1,j)
        Sim1 = Si - d1m(i,j)               ! Sim1 = tracer(i-1,j)
        Sip1 = Si + d1p(i,j)               ! Sip1 = tracer(i+1,j)

        z = d1m(i,j) - d1mm(i,j)           ! z = del2(i-1,j)
        w = d1p(i,j) - d1m(i,j)            ! w = del2(i,j)
        x = 4. * w - z
        y = 4. * z - w
        dM4m = max(0., min(x,y,z,w)) + min(0., max(x,y,z,w))
        z = d1pp(i,j) - d1p(i,j)           ! z = del2(i+1,j)
        x = 4. * w - z
        y = 4. * z - w
        dM4p = max(0., min(x,y,z,w)) + min(0., max(x,y,z,w))

        qAV = 0.5 * ( Si + Sip1 )
        x = d1p(i,j)                       ! = Sip1 - Si
        y = 3. * d1m(i,j)                  ! = 3. * ( Si - Sim1 )
        qMP = Si + max(0., min(x, y)) + min(0., max(x,y))
        qUL = Si + y
        qLC = ( Si + 0.5 * d1m(i,j) ) + fourThirds * dM4m
        qMD = qAV - 0.5 * dM4p
        qMin = max( min(qMD,Si,Sip1), min(Si,qUL,qLC) )
        qMax = min( max(qMD,Si,Sip1), max(Si,qUL,qLC) )
        if ( (aR(i,j)-Si)*(aR(i,j)-qMP).gt.1.e-10 )  &
             aR(i,j) = min( max( aR(i,j), qMin ), qMax )

        qAV = 0.5 * ( Si + Sim1 )
        x = -d1m(i,j)                      ! = Sim1 - Si
        y = -3. * d1p(i,j)                 ! = 3. * ( Si - Sip1 )
        qMP = Si + max(0., min(x, y)) + min(0., max(x,y))
        qUL = Si + y
        qLC = ( Si - 0.5 * d1p(i,j) ) + fourThirds * dM4p
        qMD = qAV - 0.5 * dM4m
        qMin = max( min(qMD,Si,Sim1), min(Si,qUL,qLC) )
        qMax = min( max(qMD,Si,Sim1), max(Si,qUL,qLC) )
        if ( (aL(i,j)-Si)*(aL(i,j)-qMP).gt.1.e-10 )  &
             aL(i,j) = min( max( aL(i,j), qMin ), qMax )
      enddo !i
    enddo !j

end subroutine ppm_limit_sh
! </SUBROUTINE> NAME="ppm_limit_sh"

! ============================================================================


subroutine xTracerFlx_ppmNiki_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_ppmNiki_cne(sv,ppm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise parabolic tracer flux in x direction.
  !             The choice of the limiter is as follows:
  !             ppm_hlimiter=0 : No limiter 

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference: Experimental. Trying to reproduce ppm with the formulation of pqm 

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  integer, intent(in) :: ppm_hlimiter           ! Choice of the limiter function
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: cfl, massflux, mcflp1
  real,dimension(INH_,JNH_,0:4)  :: poly
  real,dimension(INH_)  :: h
  real,dimension(INH_,1:2) :: edge_values,edge_slopes
  real :: uBar,uL,uR,uLprime,uRprime

  tracer_flux = 0.0
  poly = 0.0

  do j=JNLO_,JNHI_

     edge_values = 0.0 
     edge_slopes = 0.0

     !Uniform grids
     h(:)=dx(IS_ , j)
     call edge_values_explicit_h4_uniformGrid(IS_-2,IE_+2, sv(1,5), tracer_field(:,j), edge_values )


     do i=IS_-1,IE_+1

        !Note: uLprime and uRprime should contain a width, i.e., 
        ![uLprime] = [uRprime] = [tracer_field]
        !
        uBar    = tracer_field(i,j)*tmask(i,j)
        uL      = edge_values(i,1) *tmask(i,j)*tmask(i-1,j)
        uR      = edge_values(i,2) *tmask(i,j)*tmask(i+1,j)

        !White & Adcroft JCP 227 (2008)
        poly(i,j,0) = uL
        poly(i,j,1) = 6*uBar-4*uL-2*uR
        poly(i,j,2) = 3*(uL+uR-2*uBar)

     enddo !i
  enddo !j

  do j=JNLO_,JNHI_
     do I=IS_,IE_

        massflux = dy(I,j) * uhrho(I,j)

        cfl = u(I,j) * dt /dx(I,j) 
        mcflp1 = -cfl + 1.0

        if (massflux.gt.0.0) then                              
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                *(poly(i,j,0) +  poly(i,j,1)*(2.0-cfl) / 2  & !plm 
                +(poly(i,j,2)*(1.0-mcflp1**3) / 3 ) / cfl )

        else                                                  
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                *( poly(i+1,j,0) - poly(i+1,j,1)* cfl / 2    &   !plm 
                +poly(i+1,j,2)* cfl**2 / 3 )

        endif

     enddo !i                                                
  enddo !j



end subroutine xTracerFlx_ppmNiki_cne

subroutine xTracerFlx_pqm_cne(sv,pqm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_pqm_cne(sv,pqm_hlimiter,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise quartic tracer flux in x direction.
  !             The choice of the limiter is as follows:
  !             pqm_hlimiter=0 : No limiter 
  !             pqm_hlimiter=1 : 

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  integer, intent(in) :: pqm_hlimiter           ! Choice of the limiter function
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,JN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: cfl, massflux, mcflp1
  real,dimension(INH_,JNH_,0:4)  :: poly
  real,dimension(INH_)  :: h
  real,dimension(INH_,1:2) :: edge_values,edge_slopes
  real :: uBar,uL,uR,uLprime,uRprime

  tracer_flux = 0.0
  poly = 0.0

  do j=JNLO_,JNHI_

     edge_values = 0.0 
     edge_slopes = 0.0

     h(:) = dx(:,j)

     !Non-uniform grids
     !call edge_slopes_implicit_h3(IS_-2, IE_+2, sv(1,5), h, tracer_field(:,j), edge_slopes )
     !call edge_values_explicit_h4(IS_-2, IE_+2, sv(1,5), h, tracer_field(:,j), edge_values )
     !
     !nnz: NOTE: I think we have a problem here. In order for the following calls to work properly 
     !           dx should also have a halo of at least 3, which it doesn't when passed from MOM or GOLD!
     !
     !Uniform grids
     h(:)=dx(IS_ , j)
     call edge_values_explicit_h4_uniformGrid(IS_-2,IE_+2, sv(1,5), tracer_field(:,j), edge_values )
     call edge_slopes_explicit_h4_uniformGrid(IS_-2,IE_+2, sv(1,5), h(IS_), tracer_field(:,j), edge_slopes )


     ! Limit the edge values, aL and aR, for monotonicity
     if (pqm_hlimiter.eq.1) then
        call pqm_limiter(IS_,IE_+1, sv(1,5), sv(1,5), h, tracer_field(:,j), edge_slopes, edge_values)
     endif

     do I=IS_,IE_+1

        !Note: uLprime and uRprime should contain a width, i.e., 
        ![uLprime] = [uRprime] = [tracer_field]
        !
        uBar    = tracer_field(i,j)*tmask(i,j)
        uL      = edge_values(i,1) *tmask(i,j)*tmask(i-1,j)
        uLprime = edge_slopes(i,1) *tmask(i,j)*tmask(i-1,j) 
        uR      = edge_values(i,2) *tmask(i,j)*tmask(i+1,j)
        uRprime = edge_slopes(i,2) *tmask(i,j)*tmask(i+1,j)

        !White & Adcroft JCP 227 (2008)
        poly(i,j,0) = uL
        poly(i,j,1) = uLprime
        poly(i,j,2) = 30*uBar - 12*uR -18*uL + 3*(uRprime - 3*uLprime)/2 !3*(uL+uR-2*uBar)!
        poly(i,j,3) = -60*uBar + 6*uLprime -4*uRprime +28*uR + 32*uL
        poly(i,j,4) =  30*uBar +5*(uRprime - uLprime)/2 - 15*(uL+uR)

     enddo !i
  enddo !j

  do j=JNLO_,JNHI_
     do I=IS_,IE_

        massflux = dy(I,j) * uhrho(I,j)

        cfl = u(I,j) * dt /dx(I,j) 
        mcflp1 = -cfl + 1.0

        if (massflux.gt.0.0) then                              
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                *(poly(i,j,0) +  poly(i,j,1)*(2.0-cfl) / 2  & !plm 
                +(poly(i,j,2)*(1.0-mcflp1**3) / 3 &
                +poly(i,j,3)*(1.0-mcflp1**4) / 4 &
                +poly(i,j,4)*(1.0-mcflp1**5) / 5 ) / cfl )

        else                                                  
           tracer_flux(I,j) = massflux * tmask(i,j) * tmask(i+1,j) &
                *( poly(i+1,j,0) - poly(i+1,j,1)* cfl / 2    &   !plm 
                +poly(i+1,j,2)* cfl**2 / 3 &
                -poly(i+1,j,3)* cfl**3 / 4 &
                +poly(i+1,j,4)* cfl**4 / 5 )

        endif

     enddo !i                                                
  enddo !j



end subroutine xTracerFlx_pqm_cne

subroutine yTracerFlx_pqm_cne(sv,pqm_hlimiter,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_pqm_cne(sv,pqm_hlimiter,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !DESCRIPTION: Calculates the piecewise quartic tracer flux in y direction.
  !             The choice of the limiter is as follows:
  !             pqm_hlimiter=0 : No limiter 
  !             pqm_hlimiter=1 : 

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  integer, intent(in) :: pqm_hlimiter           ! Choice of the limiter function
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of face at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of cell at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)  ! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real	  :: cfl, massflux, mcflp1
  real,dimension(INH_,JNH_,0:4)  :: poly
  real,dimension(JNH_)  :: h
  real,dimension(JNH_ , 1:2) :: edge_values,edge_slopes
  real :: uBar,uL,uR,uLprime,uRprime

  tracer_flux = 0.0
  poly = 0.0

  do i=IS_,IE_

     edge_values = 0.0 
     edge_slopes = 0.0
     h(:) = dx(i,:)

     !Non-uniform grids
     !call edge_slopes_implicit_h3(JS_-2, JE_+2, sv(2,1), dy(i,:) , tracer_field(i,:), edge_slopes )
     !call edge_values_explicit_h4(JS_-2, JE_+2, sv(2,1), dy(i,:) , tracer_field(i,:), edge_values )
     !nnz: NOTE: I think we have a problem here. In order for the following calls to work properly 
     !           dx should also have a halo of at least 3, which it doesn't when passed to the current sub!
     !
     !Uniform grids
     call edge_slopes_explicit_h4_uniformGrid(JS_-2,JE_+2, sv(2,5), h(JS_), tracer_field(i,:), edge_slopes )
     call edge_values_explicit_h4_uniformGrid(JS_-2,JE_+2, sv(2,5), tracer_field(i,:), edge_values )

  ! Limit the edge values, aL and aR, for monotonicity
      if (pqm_hlimiter.eq.1) then
         call pqm_limiter(JS_,JE_, sv(2,5), sv(2,5),dy(i,:) , tracer_field(i,:), edge_slopes, edge_values)
      endif

     do J=JS_,JE_+1 
          
        uBar    = tracer_field(i,j)*tmask(i,j)
        uL      = edge_values(j,1) *tmask(i,j)*tmask(i,j-1)
        uLprime = edge_slopes(j,1) *tmask(i,j)*tmask(i,j-1)
        uR      = edge_values(j,2) *tmask(i,j)*tmask(i,j+1)
        uRprime = edge_slopes(j,2) *tmask(i,j)*tmask(i,j+1)

       !White & Adcroft JCP 227 (2008)
        poly(i,j,0) = uL
        poly(i,j,1) = uLprime
        poly(i,j,2) =  30*uBar - 12*uR -18*uL + 3*(uRprime - 3*uLprime)/2
        poly(i,j,3) = -60*uBar + 6*uLprime -4*uRprime +28*uR + 32*uL
        poly(i,j,4) =  30*uBar +5*(uRprime - uLprime)/2 - 15*(uL+uR)
        
     enddo !i
  enddo !j

  do J=JS_,JE_
     do i=IS_,IE_                                          

        massflux = dx(i,J) * vhrho(i,J)

        cfl = v(i,J) * dt /dy(i,J) 
        mcflp1 = -cfl + 1.0
 
        if (massflux.gt.0.0) then                              
          tracer_flux(i,J) = massflux * tmask(i,j) * tmask(i,j+1) &
                                *(poly(i,j,0) +  poly(i,j,1)*(2.0-cfl) / 2  & !plm 
                                +(poly(i,j,2)*(1.0-mcflp1**3) / 3 &
                                 +poly(i,j,3)*(1.0-mcflp1**4) / 4 &
                                 +poly(i,j,4)*(1.0-mcflp1**5) / 5 ) / cfl )
        else                                                  
           tracer_flux(i,J) = massflux * tmask(i,j) * tmask(i,j+1) &
                               *( poly(i,j+1,0) - poly(i,j+1,1)* cfl / 2    &   !plm 
                                 +poly(i,j+1,2)* cfl**2 / 3 &
                                 -poly(i,j+1,3)* cfl**3 / 4 &
                                 +poly(i,j+1,4)* cfl**4 / 5 )
        endif

     enddo !i                                                
  enddo !j

end subroutine yTracerFlx_pqm_cne

subroutine edge_slopes_values_plm (nsc,nec,nsh, cell_values, edge_slopes,edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_slopes(nsh:,1:), edge_values(nsh:,1:)

  integer :: i
  real :: a

  do i=nsc,nec 
     !GOLD's algorithm for calculating slopes.
     if (ABS(cell_values(i+1)-cell_values(i)) < ABS(cell_values(i)-cell_values(i-1))) then
        a= 4.0*(cell_values(i+1)-cell_values(i))
     else
        a= 4.0*(cell_values(i)-cell_values(i-1))
     endif
     if ((cell_values(i+1)-cell_values(i)) * (cell_values(i)-cell_values(i-1)) < 0.0) then
        edge_slopes(i,1)= 0.0
     elseif (ABS(cell_values(i+1)-cell_values(i-1))<ABS(a)) then
        edge_slopes(i,1)= 0.5*(cell_values(i+1)-cell_values(i-1))
     else
        edge_slopes(i,1)= 0.5*a
     endif
     edge_values(i,1) = cell_values(i) - edge_slopes(i,1)/2
  enddo
     
end subroutine edge_slopes_values_plm


subroutine edge_slopes_ppm_uniformGrid(nsc,nec,nsh, h, cell_values, edge_values, edge_slopes )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: h 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(in)  :: edge_values(nsh:,1:)
  real,    intent(out) :: edge_slopes(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge slopes based on reproducing the ppm
! -----------------------------------------------------------------------------


  do i=nsc,nec
     edge_slopes(i,1)  =  6*cell_values(i) - 4*edge_values(i,1) - 2*edge_values(i,2) 
     edge_slopes(i,2)  = -6*cell_values(i) + 2*edge_values(i,1) + 4*edge_values(i,2)  
  enddo

end subroutine edge_slopes_ppm_uniformGrid

subroutine edge_slopes_explicit_h4_uniformGrid(nsc,nec,nsh, h, cell_values, edge_slopes )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: h 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_slopes(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge slopes based on 4th-order explicit estimates on uniform grids.
! (White & Adcroft, JCP 2008 Eq. 14,37)
! -----------------------------------------------------------------------------

!  edge_slopes(nsc,1) = (-11*cell_values(nsc)+45*cell_values(nsc+1)-69*cell_values(nsc+2)+35*cell_values(nsc+3))/12/h

  do i=nsc+1,nec
!my calc     edge_slopes(i,1)  = (27*(cell_values(i)  -cell_values(i-1)) -  (cell_values(i+1)-cell_values(i-2)))/24/h
     edge_slopes(i,1) = (15*(cell_values(i)  -cell_values(i-1)) -  (cell_values(i+1)-cell_values(i-2)))/12/h
     edge_slopes(i-1,2)= edge_slopes(i,1) 
  enddo
   edge_slopes(nec-1,2)= edge_slopes(nec,1) 

!  edge_slopes(nec,1)= (-35*cell_values(nec-3)+69*cell_values(nec-2)-45*cell_values(nec+1)+11*cell_values(nec))/12/h
!  edge_slopes(nec,2)=   edge_slopes(nec+1,1) 
  
end subroutine edge_slopes_explicit_h4_uniformGrid

subroutine edge_values_explicit_h3_uniformGrid(nsc,nec,nsh, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge values based on 4th-order explicit estimates on uniform grids.
! (White & Adcroft, JCP 2008 Eq. 9)
! -----------------------------------------------------------------------------

  do i=nsc,nec
     edge_values(i,1)= ( 2*cell_values(i-1)+5*cell_values(i)-  cell_values(i+1))/6
     edge_values(i,2)= (  -cell_values(i-1)+5*cell_values(i)+2*cell_values(i+1))/6
  enddo

  
end subroutine edge_values_explicit_h3_uniformGrid

subroutine edge_values_explicit_h4_uniformGrid(nsc,nec,nsh, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge values based on 4th-order explicit estimates on uniform grids.
! (White & Adcroft, JCP 2008 Eq. 10,35,36)
! -----------------------------------------------------------------------------

!  edge_values(nsc,1)  = (25*cell_values(nsc)-23*cell_values(nsc+1)+13*cell_values(nsc+2)-3*cell_values(nsc+3))/12
!  edge_values(nsc+1,1)= ( 3*cell_values(nsc)+13*cell_values(nsc+1)- 5*cell_values(nsc+2)+  cell_values(nsc+3))/12
!  edge_values(nsc,2)  = edge_values(nsc+1,1)

  do i=nsc+1,nec
!my calc     edge_values(i,1)  = (-cell_values(i-2)+9*cell_values(i-1)+9*cell_values(i)  -cell_values(i+1))/16
    edge_values(i,1)  = (-cell_values(i-2)+7*cell_values(i-1)+7*cell_values(i)  -cell_values(i+1))/12
     edge_values(i-1,2)= edge_values(i,1)
  enddo

!  edge_values(nec,1)  = (   cell_values(nec-3)-5* cell_values(nec-2)+13*cell_values(nec-1)+3* cell_values(nec))/12
! edge_values(nec+1,1)= (-3*cell_values(nec-3)+13*cell_values(nec-2)-23*cell_values(nec-1)+25*cell_values(nec))/12
! edge_values(nec,2)  = edge_values(nec+1,1)  

  edge_values(nec-1,2)= edge_values(nec,1)
  
end subroutine edge_values_explicit_h4_uniformGrid

subroutine edge_slopes_explicit_h5_uniformGrid(nsc,nec,nsh, h, cell_values, edge_slopes )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: h 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_slopes(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge slopes based on 5th-order explicit estimates on uniform grids.
! (White & Adcroft, JCP 2008 Eq. 15)
! -----------------------------------------------------------------------------

  do i=nsc+2,nec-2
     edge_slopes(i,1)= (245*(cell_values(i)  -cell_values(i-1)) - 25*(cell_values(i+1)-cell_values(i-2))&
                        +2* (cell_values(i+2)-cell_values(i-3)))/180/h
     edge_slopes(i,2)= (245*(cell_values(i+1)-cell_values(i)  ) - 25*(cell_values(i+2)-cell_values(i-1))&
                        +2* (cell_values(i+3)-cell_values(i-2)))/180/h
  enddo
  
end subroutine edge_slopes_explicit_h5_uniformGrid
 
subroutine edge_values_explicit_h5_uniformGrid(nsc,nec,nsh, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
  !
  integer :: i
! -----------------------------------------------------------------------------
! Compute edge values based on 5th-order explicit estimates on uniform grids.
! (White & Adcroft, JCP 2008 Eq. 11)
! -----------------------------------------------------------------------------

  do i=nsc+1,nec-1
     edge_values(i,1)= ( -3*cell_values(i-2)+27*cell_values(i-1)+47*cell_values(i) &
                        -13*cell_values(i+1)+2* cell_values(i+2))/60
     edge_values(i,2)= (  3*cell_values(i-2)-13*cell_values(i-1)+47*cell_values(i) &
                        +27*cell_values(i+1)-3* cell_values(i+2))/60 
  enddo
  
end subroutine edge_values_explicit_h5_uniformGrid


subroutine edge_slopes_implicit_h3(nsc,nec,nsh, h, cell_values, edge_slopes )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: h(nsh:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_slopes(nsh:,1:)
! -----------------------------------------------------------------------------
! Compute edge slopes based on third-order implicit estimates. Note that
! the estimates are fourth-order accurate on uniform grids
!
! Third-order implicit estimates of edge slopes are based on a two-cell 
! stencil. A tridiagonal system is set up and is based on expressing the 
! edge slopes in terms of neighboring cell averages. The generic 
! relationship is
!
! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} = 
! a \bar{u}_i + b \bar{u}_{i+1}
!
! and the stencil looks like this
!
!     	   i	 i+1
!	..--o------o------o--..
!     i-1/2  i+1/2  i+3/2
!
! In this routine, the coefficients \alpha, \beta, a and b are computed, 
! the tridiagonal system is built, boundary conditions are prescribed and 
! the system is solved to yield edge-slope estimates. 
! 
! There are N+1 unknowns and we are able to write N-1 equations. The 
! boundary conditions close the system.
! -----------------------------------------------------------------------------

  ! Local variables
  integer				:: i, j;				! loop indexes
  integer				:: N,offset;					! number of cells
  real					:: h0, h1;				! cell widths
  real					:: h0_2, h1_2, h0h1;
  real					:: h0_3, h1_3;
  real					:: d;
  real					:: alpha, beta;			! stencil coefficients
  real					:: a, b;
  real, dimension(5)	:: x;					! system used to enforce
  real, dimension(4,4)	:: Asys;				! boundary conditions
  real, dimension(4)	:: Bsys, Csys;
  real, dimension(3)	:: Dsys;
  real					:: eps;			
  real :: evaluation_polynomial
  real, dimension(:), allocatable :: tri_l;	! trid. system (lower diagonal)
  real, dimension(:), allocatable :: tri_d;	! trid. system (middle diagonal)
  real, dimension(:), allocatable :: tri_u;	! trid. system (upper diagonal)
  real, dimension(:), allocatable :: tri_x;	! trid. system (unknowns vector)
  real, dimension(:), allocatable :: tri_b;	! trid. system (rhs)

  ! Get number of cells (there are N+1 edge values to estimate)
  eps = 0.0 ! 1e-3;
  N = nec-nsc+1
  offset = nsc-1

  allocate ( tri_l(N+1) )	
  allocate ( tri_d(N+1) )	
  allocate ( tri_u(N+1) )	
  allocate ( tri_b(N+1) )	
  allocate ( tri_x(N+1) )	

  ! Loop on cells (except last one)
  do i = 1,N-1

	! Get cell widths
	h0 = h(offset+i) + eps;
	h1 = h(offset+i+1) + eps;

	! Auxiliary calculations
	h0h1 = h0 * h1;
	h0_2 = h0 * h0;
	h1_2 = h1 * h1;
	h0_3 = h0_2 * h0;
	h1_3 = h1_2 * h1;

	d = 4.0 * h0h1 * ( h0 + h1 ) + h1_3 + h0_3;

	! Coefficients
	alpha = h1 * (h0_2 + h0h1 - h1_2) / d;
	beta  = h0 * (h1_2 + h0h1 - h0_2) / d;
	a = -12.0 * h0h1 / d;
	b = -a;
	
	tri_l(i+1) = alpha;
	tri_d(i+1) = 1.0;
	tri_u(i+1) = beta;
	
	tri_b(i+1) = a * cell_values(offset+i) + b * cell_values(offset+i+1);
 	
  end do ! end loop on cells


  ! Boundary conditions: left boundary
  x(1) = 0.0;
  do i = 2,5
  	x(i) = x(i-1) + h(offset+i-1) + eps;
  end do

  do i = 1,4
    
	do j = 1,4
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j;
	end do
	
	Bsys(i) = cell_values(offset+i) * ( h(offset+i) + eps );
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 4 );		
  
  Dsys(1) = Csys(2); 
  Dsys(2) = 2.0 * Csys(3); 
  Dsys(3) = 3.0 * Csys(4); 
  
  tri_d(1) = 1.0;
  tri_u(1) = 0.0;
  tri_b(1) = evaluation_polynomial ( Dsys, 3, x(1) );		! first edge slope
  
  ! Boundary conditions: right boundary
  x(1) = 0.0;
  do i = 2,5
  	x(i) = x(i-1) + h(offset+N-5+i) + eps;
  end do

  do i = 1,4
    
	do j = 1,4
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j;
	end do
	
	Bsys(i) = cell_values(offset+N-4+i) * ( h(offset+N-4+i) + eps );
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 4 );		
  
  Dsys(1) = Csys(2); 
  Dsys(2) = 2.0 * Csys(3); 
  Dsys(3) = 3.0 * Csys(4); 
  
  tri_l(N+1) = 0.0;
  tri_d(N+1) = 1.0;
  tri_b(N+1) = evaluation_polynomial ( Dsys, 3, x(5) );		! last edge slope

  ! Solve tridiagonal system and assign edge values
  call solve_tridiagonal_system ( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 );

  do i = 2,N
     edge_slopes(offset+i,1)   = tri_x(i);
     edge_slopes(offset+i-1,2) = tri_x(i);
  end do
  edge_slopes(offset+1,1) = tri_x(1);
  edge_slopes(offset+N,2) = tri_x(N+1);



  deallocate(tri_l,tri_d,tri_u,tri_b,tri_x) 

!  do i=nsc,nec
!     edge_slopes(i,1)= (cell_values(i)-cell_values(i-1)) !this is a crude estimate for slope*width
!     edge_slopes(i-1,2)= edge_slopes(i,1)
!  enddo

end subroutine edge_slopes_implicit_h3


subroutine edge_slopes_implicit_h5(nsc,nec,nsd,nsh, h, cell_values, edge_slopes )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsd,nsh
  real,    intent(in)  :: h(nsd:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_slopes(nsh:,1:)
! -----------------------------------------------------------------------------
! Fifth-order implicit estimates of edge values are based on a four-cell, 
! three-edge stencil. A tridiagonal system is set up and is based on 
! expressing the edge slopes in terms of neighboring cell averages. 
! 
! The generic relationship is
!
! \alpha u'_{i-1/2} + u'_{i+1/2} + \beta u'_{i+3/2} = 
! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2}
!
! and the stencil looks like this
!
!     	  i-1     i	    i+1    i+2
!	..--o------o------o------o------o--..
!      		 i-1/2  i+1/2  i+3/2
!
! In this routine, the coefficients \alpha, \beta, a, b, c and d are 
! computed, the tridiagonal system is built, boundary conditions are 
! prescribed and the system is solved to yield edge-value estimates. 
!
! Note that the centered stencil only applies to edges 3 to N-1 (edges are
! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other
! equations are written by using a right-biased stencil for edge 2 and a
! left-biased stencil for edge N. The prescription of boundary conditions
! (using sixth-order polynomials) closes the system.
!
! CAUTION: For each edge, in order to determine the coefficients of the
! 		   implicit expression, a 6x6 linear system is solved. This may
!          become computationally expensive if regridding is carried out
!  		   often. Figuring out closed-form expressions for these coefficients
!          on nonuniform meshes turned out to be intractable.
! -----------------------------------------------------------------------------

  ! Local variables
  integer				:: i, j, k 				! loop indexes
  integer				:: N, offset 					! number of cells
  real					:: h0, h1, h2, h3 		! cell widths
  real					:: g, g_2, g_3 			! the following are	
  real					:: g_4, g_5, g_6 		! auxiliary variables
  real					:: d2, d3, d4, d5, d6   ! to set up the systems
  real					:: n2, n3, n4, n5, n6 	! used to compute the
  real					:: h1_2, h2_2 			! the coefficients of the
  real					:: h1_3, h2_3 			! tridiagonal system
  real					:: h1_4, h2_4 			! ...
  real					:: h1_5, h2_5 			! ...
  real					:: h1_6, h2_6 			! ...
  real					:: h0ph1, h0ph1_2 		! ...
  real					:: h0ph1_3, h0ph1_4 	! ...
  real					:: h2ph3, h2ph3_2 		! ...
  real					:: h2ph3_3, h2ph3_4 	! ...
  real					:: alpha, beta 			! stencil coefficients
  real					:: a, b, c, d 			! "
  real, dimension(7)	:: x 					! system used to enforce
  real, dimension(6,6)	:: Asys 				! boundary conditions
  real, dimension(6)	:: Bsys, Csys 			! ...
  real, dimension(5)	:: Dsys 				! derivative		
  real					:: eps 					! 
  real :: evaluation_polynomial
  real, dimension(:), allocatable :: tri_l;	! trid. system (lower diagonal)
  real, dimension(:), allocatable :: tri_d;	! trid. system (middle diagonal)
  real, dimension(:), allocatable :: tri_u;	! trid. system (upper diagonal)
  real, dimension(:), allocatable :: tri_x;	! trid. system (unknowns vector)
  real, dimension(:), allocatable :: tri_b;	! trid. system (rhs)

  ! Get number of cells (there are N+1 edge values to estimate)
  eps = 0.0 !1e-2 

  N = nec-nsc+1
  offset = nsc-1

  allocate ( tri_l(N+1) )	
  allocate ( tri_d(N+1) )	
  allocate ( tri_u(N+1) )	
  allocate ( tri_b(N+1) )	
  allocate ( tri_x(N+1) )	

  ! Loop on cells (except last one)
  do k = 2,N-2

	! Cell widths
    h0 = h(offset+k-1) + eps 
    h1 = h(offset+k+0) + eps 
    h2 = h(offset+k+1) + eps 
    h3 = h(offset+k+2) + eps 
	
	! Auxiliary calculations
	h1_2 = h1 * h1 
	h1_3 = h1_2 * h1  
	h1_4 = h1_2 * h1_2 
	h1_5 = h1_3 * h1_2 
	h1_6 = h1_3 * h1_3 
	
	h2_2 = h2 * h2 
	h2_3 = h2_2 * h2  
	h2_4 = h2_2 * h2_2 
	h2_5 = h2_3 * h2_2 
	h2_6 = h2_3 * h2_3 
	
	g   = h0 + h1 
	g_2 = g * g 
	g_3 = g * g_2 
	g_4 = g_2 * g_2 
	g_5 = g_4 * g 
	g_6 = g_3 * g_3 

	d2 = ( h1_2 - g_2 ) / h0 
	d3 = ( h1_3 - g_3 ) / h0 
	d4 = ( h1_4 - g_4 ) / h0 
	d5 = ( h1_5 - g_5 ) / h0 
	d6 = ( h1_6 - g_6 ) / h0 
	
	g   = h2 + h3 
	g_2 = g * g 
	g_3 = g * g_2 
	g_4 = g_2 * g_2 
	g_5 = g_4 * g 
	g_6 = g_3 * g_3 
	
	n2 = ( g_2 - h2_2 ) / h3 
	n3 = ( g_3 - h2_3 ) / h3 
	n4 = ( g_4 - h2_4 ) / h3 
	n5 = ( g_5 - h2_5 ) / h3 
	n6 = ( g_6 - h2_6 ) / h3 

	! Compute matrix entries
	Asys(1,1) =	0.0 
	Asys(1,2) =	0.0 
	Asys(1,3) =	1.0 
	Asys(1,4) =	1.0 
	Asys(1,5) =	1.0 
	Asys(1,6) =	1.0 
	
	Asys(2,1) = 1.0 
	Asys(2,2) = 1.0 
	Asys(2,3) =	-0.5 * d2 
	Asys(2,4) =	0.5 * h1 
	Asys(2,5) =	-0.5 * h2 
	Asys(2,6) =	-0.5 * n2 
	
	Asys(3,1) =	h1 
	Asys(3,2) = - h2 
	Asys(3,3) =	- d3 / 6.0 
	Asys(3,4) =	h1_2 / 6.0 
	Asys(3,5) =	h2_2 / 6.0 
	Asys(3,6) =	n3 / 6.0 
	
	Asys(4,1) =	- h1_2 / 2.0 
	Asys(4,2) = - h2_2 / 2.0 
	Asys(4,3) =	d4 / 24.0 
	Asys(4,4) =	- h1_3 / 24.0 
	Asys(4,5) =	h2_3 / 24.0 
	Asys(4,6) =	n4 / 24.0 
	
	Asys(5,1) =	h1_3 / 6.0 
	Asys(5,2) = - h2_3 / 6.0 
	Asys(5,3) =	- d5 / 120.0 
	Asys(5,4) =	h1_4 / 120.0 
	Asys(5,5) =	h2_4 / 120.0 
	Asys(5,6) =	n5 / 120.0 
	
	Asys(6,1) =	- h1_4 / 24.0 
	Asys(6,2) = - h2_4 / 24.0 
	Asys(6,3) =	d6 / 720.0 
	Asys(6,4) =	- h1_5 / 720.0 
	Asys(6,5) =	h2_5 / 720.0 
	Asys(6,6) =	n6 / 720.0 
	
	Bsys(:) = (/ 0.0, -1.0, 0.0, 0.0, 0.0, 0.0 /) 
 
    call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

	alpha = Csys(1) 
	beta  = Csys(2) 
	a = Csys(3) 
	b = Csys(4) 
	c = Csys(5) 
	d = Csys(6) 

	tri_l(k+1) = alpha 
	tri_d(k+1) = 1.0 
	tri_u(k+1) = beta 
	tri_b(k+1) = a * cell_values(offset+k-1) + b * cell_values(offset+k) + c * cell_values(offset+k+1) + d * cell_values(offset+k+2) 
 	
  end do ! end loop on cells

  ! Use a right-biased stencil for the second row
  
  ! Cell widths
  h0 = h(offset+1) + eps 
  h1 = h(offset+2) + eps 
  h2 = h(offset+3) + eps 
  h3 = h(offset+4) + eps 
	
  ! Auxiliary calculations
  h1_2 = h1 * h1 
  h1_3 = h1_2 * h1  
  h1_4 = h1_2 * h1_2 
  h1_5 = h1_3 * h1_2 
  h1_6 = h1_3 * h1_3 
	
  h2_2 = h2 * h2 
  h2_3 = h2_2 * h2  
  h2_4 = h2_2 * h2_2 
  h2_5 = h2_3 * h2_2 
  h2_6 = h2_3 * h2_3 
	
  g   = h0 + h1 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 

  h0ph1   = h0 + h1 
  h0ph1_2 = h0ph1 * h0ph1  
  h0ph1_3 = h0ph1_2 * h0ph1  
  h0ph1_4 = h0ph1_2 * h0ph1_2  

  d2 = ( h1_2 - g_2 ) / h0 
  d3 = ( h1_3 - g_3 ) / h0 
  d4 = ( h1_4 - g_4 ) / h0 
  d5 = ( h1_5 - g_5 ) / h0 
  d6 = ( h1_6 - g_6 ) / h0 
	
  g   = h2 + h3 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 
	
  n2 = ( g_2 - h2_2 ) / h3 
  n3 = ( g_3 - h2_3 ) / h3 
  n4 = ( g_4 - h2_4 ) / h3 
  n5 = ( g_5 - h2_5 ) / h3 
  n6 = ( g_6 - h2_6 ) / h3 

  ! Compute matrix entries
  Asys(1,1) = 0.0 
  Asys(1,2) = 0.0 
  Asys(1,3) = 1.0 
  Asys(1,4) = 1.0 
  Asys(1,5) = 1.0 
  Asys(1,6) = 1.0 
	
  Asys(2,1) = 1.0 
  Asys(2,2) = 1.0 
  Asys(2,3) = -0.5 * d2 
  Asys(2,4) = 0.5 * h1 
  Asys(2,5) = -0.5 * h2 
  Asys(2,6) = -0.5 * n2 
	
  Asys(3,1) = h0ph1 
  Asys(3,2) = 0.0 
  Asys(3,3) = - d3 / 6.0 
  Asys(3,4) = h1_2 / 6.0 
  Asys(3,5) = h2_2 / 6.0 
  Asys(3,6) = n3 / 6.0 

  Asys(4,1) = - h0ph1_2 / 2.0 
  Asys(4,2) = 0.0 
  Asys(4,3) = d4 / 24.0 
  Asys(4,4) = - h1_3 / 24.0 
  Asys(4,5) = h2_3 / 24.0 
  Asys(4,6) = n4 / 24.0 

  Asys(5,1) = h0ph1_3 / 6.0 
  Asys(5,2) = 0.0 
  Asys(5,3) = - d5 / 120.0 
  Asys(5,4) = h1_4 / 120.0 
  Asys(5,5) = h2_4 / 120.0 
  Asys(5,6) = n5 / 120.0 

  Asys(6,1) = - h0ph1_4 / 24.0 
  Asys(6,2) = 0.0 
  Asys(6,3) = d6 / 720.0 
  Asys(6,4) = - h1_5 / 720.0 
  Asys(6,5) = h2_5 / 720.0 
  Asys(6,6) = n6 / 720.0 

  Bsys(:) = (/ 0.0, -1.0, -h1, h1_2/2.0, -h1_3/6.0, h1_4/24.0 /) 

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

  alpha = Csys(1) 
  beta  = Csys(2) 
  a = Csys(3) 
  b = Csys(4) 
  c = Csys(5) 
  d = Csys(6) 

  tri_l(2) = alpha 
  tri_d(2) = 1.0 
  tri_u(2) = beta 
  tri_b(2) = a * cell_values(offset+1) + b * cell_values(offset+2) + c * cell_values(offset+3) + d * cell_values(offset+4) 
  
  ! Boundary conditions: left boundary
  x(1) = 0.0 
  do i = 2,7
  	x(i) = x(i-1) + h(offset+i-1) + eps 
  end do

  do i = 1,6
    
	do j = 1,6
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j 
	end do
	
	Bsys(i) = cell_values(offset+i) * ( h(offset+i) + eps ) 
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		
  
  Dsys(1) = Csys(2)  
  Dsys(2) = 2.0 * Csys(3)  
  Dsys(3) = 3.0 * Csys(4)  
  Dsys(4) = 4.0 * Csys(5)  
  Dsys(5) = 5.0 * Csys(6)  
  
  tri_d(1) = 0.0 
  tri_d(1) = 1.0 
  tri_u(1) = 0.0 
  tri_b(1) = evaluation_polynomial ( Dsys, 5, x(1) ) 		! first edge value
  
  ! Use a left-biased stencil for the second to last row
  
  ! Cell widths
  h0 = h(offset+N-3) + eps 
  h1 = h(offset+N-2) + eps 
  h2 = h(offset+N-1) + eps 
  h3 = h(offset+N) + eps 
	
  ! Auxiliary calculations
  h1_2 = h1 * h1 
  h1_3 = h1_2 * h1  
  h1_4 = h1_2 * h1_2 
  h1_5 = h1_3 * h1_2 
  h1_6 = h1_3 * h1_3 
	
  h2_2 = h2 * h2 
  h2_3 = h2_2 * h2  
  h2_4 = h2_2 * h2_2 
  h2_5 = h2_3 * h2_2 
  h2_6 = h2_3 * h2_3 
	
  g   = h0 + h1 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 

  h2ph3   = h2 + h3 
  h2ph3_2 = h2ph3 * h2ph3  
  h2ph3_3 = h2ph3_2 * h2ph3  
  h2ph3_4 = h2ph3_2 * h2ph3_2  

  d2 = ( h1_2 - g_2 ) / h0 
  d3 = ( h1_3 - g_3 ) / h0 
  d4 = ( h1_4 - g_4 ) / h0 
  d5 = ( h1_5 - g_5 ) / h0 
  d6 = ( h1_6 - g_6 ) / h0 
	
  g   = h2 + h3 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 
	
  n2 = ( g_2 - h2_2 ) / h3 
  n3 = ( g_3 - h2_3 ) / h3 
  n4 = ( g_4 - h2_4 ) / h3 
  n5 = ( g_5 - h2_5 ) / h3 
  n6 = ( g_6 - h2_6 ) / h3 

  ! Compute matrix entries
  Asys(1,1) = 0.0 
  Asys(1,2) = 0.0 
  Asys(1,3) = 1.0 
  Asys(1,4) = 1.0 
  Asys(1,5) = 1.0 
  Asys(1,6) = 1.0 
	
  Asys(2,1) = 1.0 
  Asys(2,2) = 1.0 
  Asys(2,3) = -0.5 * d2 
  Asys(2,4) = 0.5 * h1 
  Asys(2,5) = -0.5 * h2 
  Asys(2,6) = -0.5 * n2 
	
  Asys(3,1) = 0.0 
  Asys(3,2) = - h2ph3 
  Asys(3,3) = - d3 / 6.0 
  Asys(3,4) = h1_2 / 6.0 
  Asys(3,5) = h2_2 / 6.0 
  Asys(3,6) = n3 / 6.0 

  Asys(4,1) = 0.0 
  Asys(4,2) = - h2ph3_2 / 2.0 
  Asys(4,3) = d4 / 24.0 
  Asys(4,4) = - h1_3 / 24.0 
  Asys(4,5) = h2_3 / 24.0 
  Asys(4,6) = n4 / 24.0 

  Asys(5,1) = 0.0 
  Asys(5,2) = - h2ph3_3 / 6.0 
  Asys(5,3) = - d5 / 120.0 
  Asys(5,4) = h1_4 / 120.0 
  Asys(5,5) = h2_4 / 120.0 
  Asys(5,6) = n5 / 120.0 

  Asys(6,1) = 0.0 
  Asys(6,2) = - h2ph3_4 / 24.0 
  Asys(6,3) = d6 / 720.0 
  Asys(6,4) = - h1_5 / 720.0 
  Asys(6,5) = h2_5 / 720.0 
  Asys(6,6) = n6 / 720.0 

  Bsys(:) = (/ 0.0, -1.0, h2, h2_2/2.0, h2_3/6.0, h2_4/24.0 /) 

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

  alpha = Csys(1) 
  beta  = Csys(2) 
  a = Csys(3) 
  b = Csys(4) 
  c = Csys(5) 
  d = Csys(6) 

  tri_l(N) = alpha 
  tri_d(N) = 1.0 
  tri_u(N) = beta 
  tri_b(N) = a * cell_values(offset+N-3) + b * cell_values(offset+N-2) + c * cell_values(offset+N-1) + d * cell_values(offset+N) 
  
  ! Boundary conditions: right boundary
  x(1) = 0.0 
  do i = 2,7
  	x(i) = x(i-1) + h(offset+N-7+i) + eps 
  end do

  do i = 1,6
    
	do j = 1,6
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j 
	end do
	
	Bsys(i) = cell_values(offset+N-6+i) * ( h(offset+N-6+i) + eps ) 
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		
  
  Dsys(1) = Csys(2)  
  Dsys(2) = 2.0 * Csys(3)  
  Dsys(3) = 3.0 * Csys(4)  
  Dsys(4) = 4.0 * Csys(5)  
  Dsys(5) = 5.0 * Csys(6)  
  
  tri_l(N+1) = 0.0 
  tri_d(N+1) = 1.0 
  tri_u(N+1) = 0.0 
  tri_b(N+1) = evaluation_polynomial ( Dsys, 5, x(7) ) 		! last edge value

  ! Solve tridiagonal system and assign edge values
  call solve_tridiagonal_system ( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) 

  do i = 2,N
     edge_slopes(offset+i,1)   = tri_x(i) 
     edge_slopes(offset+i-1,2) = tri_x(i) 
  end do
  edge_slopes(offset+1,1) = tri_x(1) 
  edge_slopes(offset+N,2) = tri_x(N+1) 


  deallocate(tri_l,tri_d,tri_u,tri_b,tri_x) 


end subroutine edge_slopes_implicit_h5


!------------------------------------------------------------------------------
! Compute h4 edge values (explicit fourth order accurate)
!------------------------------------------------------------------------------
subroutine edge_values_explicit_h4 (nsc,nec,nsh, h, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: h(nsh:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
! -----------------------------------------------------------------------------
! Compute edge values based on fourth-order explicit estimates.
! These estimates are based on a cubic interpolant spanning four cells 
! and evaluated at the location of the middle edge. An interpolant spanning 
! cells i-2, i-1, i and i+1 is evaluated at edge i-1/2. The estimate for 
! each edge is unique.
!
!       i-2    i-1     i     i+1
! ..--o------o------o------o------o--..
!				  i-1/2
! 
! The first two edge values are estimated by evaluating the first available 
! cubic interpolant, i.e., the interpolant spanning cells 1, 2, 3 and 4. 
! Similarly, the last two edge values are estimated by evaluating the last 
! available interpolant.
!
! For this fourth-order scheme, at least four cells must exist.
! -----------------------------------------------------------------------------

  ! Local variables
  integer				:: i, j;
  integer				:: N,offset;				! number of grid cells
  real					:: u0, u1, u2, u3;
  real					:: h0, h1, h2, h3;
  real					:: f1, f2, f3;		! auxiliary variables
  real					:: e;				! edge value	
  real, dimension(5)	:: x;				! used to compute edge
  real, dimension(4,4)	:: A;				! values near the boundaries
  real, dimension(4)	:: B, C;
  real					:: eps;			
  real :: evaluation_polynomial
  real, dimension(:), allocatable :: tri_l;	! trid. system (lower diagonal)
  real, dimension(:), allocatable :: tri_d;	! trid. system (middle diagonal)
  real, dimension(:), allocatable :: tri_u;	! trid. system (upper diagonal)
  real, dimension(:), allocatable :: tri_x;	! trid. system (unknowns vector)
  real, dimension(:), allocatable :: tri_b;	! trid. system (rhs)
  
  eps = 0.0! 1e-3;

  N = nec-nsc+1
  offset = nsc-1

  allocate ( tri_l(N+1) )	
  allocate ( tri_d(N+1) )	
  allocate ( tri_u(N+1) )	
  allocate ( tri_b(N+1) )	
  allocate ( tri_x(N+1) )	
 
  ! Loop on interior cells
  do i = 3,N-1
	
    h0 = h(offset+i-2) + eps;
    h1 = h(offset+i-1) + eps;
    h2 = h(offset+i) + eps;
    h3 = h(offset+i+1) + eps;
    
	u0 = cell_values(offset+i-2);
    u1 = cell_values(offset+i-1);
    u2 = cell_values(offset+i);
    u3 = cell_values(offset+i+1);
	
	f1 = (h0+h1) * (h2+h3) / (h1+h2);
	f2 = u1 * h2 + u2 * h1;
	f3 = 1.0 / (h0+h1+h2) + 1.0 / (h1+h2+h3);
	
    e = f1 * f2 * f3;      

	f1 = h2 * (h2+h3) / ( (h0+h1+h2)*(h0+h1) );
	f2 = u1*(h0+2.0*h1) - u0*h1;

	e = e + f1*f2;
	
	f1 = h1 * (h0+h1) / ( (h1+h2+h3)*(h2+h3) );
	f2 = u2*(2.0*h2+h3) - u3*h2;
	
	e = e + f1*f2;
	
	e = e / ( h0 + h1 + h2 + h3);
	
	edge_values(offset+i,1) = e;
	edge_values(offset+i-1,2) = e;

  end do ! end loop on interior cells


  ! Determine first two edge values
  x(1) = 0.0;
  do i = 2,5
  	x(i) = x(i-1) + h(offset+i-1) + eps;
  end do

  do i = 1,4
    
	do j = 1,4
	  A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j;
	end do
	
	B(i) = cell_values(offset+i) * ( h(offset+i) + eps );
	
  end do	

  call solve_linear_system ( A, B, C, 4 );		

  ! First edge value
  edge_values(offset+1,1) = evaluation_polynomial ( C, 4, x(1) );
  
  ! Second edge value
  edge_values(offset+1,2) = evaluation_polynomial ( C, 4, x(2) );
  edge_values(offset+2,1) = edge_values(offset+1,2);								

  ! Determine last two edge values
  x(1) = 0.0;
  do i = 2,5
  	x(i) = x(i-1) + h(offset+N-5+i) + eps;
  end do

  do i = 1,4
    
	do j = 1,4
	  A(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j;
	end do
	
	B(i) = cell_values(offset+N-4+i) * ( h(offset+N-4+i) + eps );
	
  end do	

  call solve_linear_system ( A, B, C, 4 );		
  
  ! Last edge value
  edge_values(offset+N,2) = evaluation_polynomial ( C, 4, x(5) );
  
  ! Second to last edge value
  edge_values(offset+N,1) = evaluation_polynomial ( C, 4, x(4) );	
  
  edge_values(offset+N-1,2) = edge_values(offset+N,1);

  deallocate(tri_l,tri_d,tri_u,tri_b,tri_x) 

!  do i=nsc,nec
!     edge_values(i,1)= (cell_values(i)+cell_values(i-1))/2
!     edge_values(i-1,2)= edge_values(i,1)
!  enddo
end subroutine edge_values_explicit_h4


!------------------------------------------------------------------------------
! Compute ih6 edge values (implicit sixth order accurate)
!------------------------------------------------------------------------------
subroutine edge_values_implicit_h6 (nsc,nec,nsd,nsh, h, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsd,nsh
  real,    intent(in)  :: h(nsd:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
! -----------------------------------------------------------------------------
! Sixth-order implicit estimates of edge values are based on a four-cell, 
! three-edge stencil. A tridiagonal system is set up and is based on 
! expressing the edge values in terms of neighboring cell averages. 
! 
! The generic relationship is
!
! \alpha u_{i-1/2} + u_{i+1/2} + \beta u_{i+3/2} = 
! a \bar{u}_{i-1} + b \bar{u}_i + c \bar{u}_{i+1} + d \bar{u}_{i+2}
!
! and the stencil looks like this
!
!     	  i-1     i	    i+1    i+2
!	..--o------o------o------o------o--..
!      		 i-1/2  i+1/2  i+3/2
!
! In this routine, the coefficients \alpha, \beta, a, b, c and d are 
! computed, the tridiagonal system is built, boundary conditions are 
! prescribed and the system is solved to yield edge-value estimates. 
!
! Note that the centered stencil only applies to edges 3 to N-1 (edges are
! numbered 1 to n+1), which yields N-3 equations for N+1 unknowns. Two other
! equations are written by using a right-biased stencil for edge 2 and a
! left-biased stencil for edge N. The prescription of boundary conditions
! (using sixth-order polynomials) closes the system.
!
! CAUTION: For each edge, in order to determine the coefficients of the
! 		   implicit expression, a 6x6 linear system is solved. This may
!          become computationally expensive if regridding is carried out
!  		   often. Figuring out closed-form expressions for these coefficients
!          on nonuniform meshes turned out to be intractable.
! -----------------------------------------------------------------------------

  ! Local variables
  integer				:: i, j, k 				! loop indexes
  integer				:: N,offset 					! number of cells
  real					:: h0, h1, h2, h3 		! cell widths
  real					:: g, g_2, g_3 			! the following are	
  real					:: g_4, g_5, g_6 		! auxiliary variables
  real					:: d2, d3, d4, d5, d6   ! to set up the systems
  real					:: n2, n3, n4, n5, n6 	! used to compute the
  real					:: h1_2, h2_2 			! the coefficients of the
  real					:: h1_3, h2_3 			! tridiagonal system
  real					:: h1_4, h2_4 			! ...
  real					:: h1_5, h2_5 			! ...
  real					:: h1_6, h2_6 			! ...
  real					:: h0ph1, h0ph1_2 		! ...
  real					:: h0ph1_3, h0ph1_4 	! ...
  real					:: h2ph3, h2ph3_2 		! ...
  real					:: h2ph3_3, h2ph3_4 	! ...
  real					:: h0ph1_5, h2ph3_5 	! ...
  real					:: alpha, beta 			! stencil coefficients
  real					:: a, b, c, d 			! "
  real, dimension(7)	:: x 					! system used to enforce
  real, dimension(6,6)	:: Asys 				! boundary conditions
  real, dimension(6)	:: Bsys, Csys 			! ...
  real					:: eps 					! 
  real :: evaluation_polynomial
  real, dimension(:), allocatable :: tri_l;	! trid. system (lower diagonal)
  real, dimension(:), allocatable :: tri_d;	! trid. system (middle diagonal)
  real, dimension(:), allocatable :: tri_u;	! trid. system (upper diagonal)
  real, dimension(:), allocatable :: tri_x;	! trid. system (unknowns vector)
  real, dimension(:), allocatable :: tri_b;	! trid. system (rhs)

  ! Get number of cells (there are N+1 edge values to estimate)
  eps = 0.0 !1e-2 


  N = nec-nsc+1
  offset = nsc-1

  allocate ( tri_l(N+1) )	
  allocate ( tri_d(N+1) )	
  allocate ( tri_u(N+1) )	
  allocate ( tri_b(N+1) )	
  allocate ( tri_x(N+1) )	

  ! Loop on cells (except last one)
  do k = 2,N-2

	! Cell widths
    h0 = h(offset+k-1) + eps 
    h1 = h(offset+k+0) + eps 
    h2 = h(offset+k+1) + eps 
    h3 = h(offset+k+2) + eps 

	! Auxiliary calculations
	h1_2 = h1 * h1 
	h1_3 = h1_2 * h1  
	h1_4 = h1_2 * h1_2 
	h1_5 = h1_3 * h1_2 
	h1_6 = h1_3 * h1_3 
	
	h2_2 = h2 * h2 
	h2_3 = h2_2 * h2  
	h2_4 = h2_2 * h2_2 
	h2_5 = h2_3 * h2_2 
	h2_6 = h2_3 * h2_3 
	
	g   = h0 + h1 
	g_2 = g * g 
	g_3 = g * g_2 
	g_4 = g_2 * g_2 
	g_5 = g_4 * g 
	g_6 = g_3 * g_3 

	d2 = ( h1_2 - g_2 ) / h0 
	d3 = ( h1_3 - g_3 ) / h0 
	d4 = ( h1_4 - g_4 ) / h0 
	d5 = ( h1_5 - g_5 ) / h0 
	d6 = ( h1_6 - g_6 ) / h0 
	
	g   = h2 + h3 
	g_2 = g * g 
	g_3 = g * g_2 
	g_4 = g_2 * g_2 
	g_5 = g_4 * g 
	g_6 = g_3 * g_3 
	
	n2 = ( g_2 - h2_2 ) / h3 
	n3 = ( g_3 - h2_3 ) / h3 
	n4 = ( g_4 - h2_4 ) / h3 
	n5 = ( g_5 - h2_5 ) / h3 
	n6 = ( g_6 - h2_6 ) / h3 

	! Compute matrix entries
	Asys(1,1) =	1.0 
	Asys(1,2) =	1.0 
	Asys(1,3) =	-1.0 
	Asys(1,4) =	-1.0 
	Asys(1,5) =	-1.0 
	Asys(1,6) =	-1.0 
	
	Asys(2,1) =	- h1 
	Asys(2,2) = h2 
	Asys(2,3) =	-0.5 * d2 
	Asys(2,4) =	0.5 * h1 
	Asys(2,5) =	-0.5 * h2 
	Asys(2,6) =	-0.5 * n2 
	
	Asys(3,1) =	0.5 * h1_2 
	Asys(3,2) = 0.5 * h2_2 
	Asys(3,3) =	d3 / 6.0 
	Asys(3,4) =	- h1_2 / 6.0 
	Asys(3,5) =	- h2_2 / 6.0 
	Asys(3,6) =	- n3 / 6.0 
	
	Asys(4,1) =	- h1_3 / 6.0 
	Asys(4,2) = h2_3 / 6.0 
	Asys(4,3) =	- d4 / 24.0 
	Asys(4,4) =	h1_3 / 24.0 
	Asys(4,5) =	- h2_3 / 24.0 
	Asys(4,6) =	- n4 / 24.0 
	
	Asys(5,1) =	h1_4 / 24.0 
	Asys(5,2) = h2_4 / 24.0 
	Asys(5,3) =	d5 / 120.0 
	Asys(5,4) =	- h1_4 / 120.0 
	Asys(5,5) =	- h2_4 / 120.0 
	Asys(5,6) =	- n5 / 120.0 
	
	Asys(6,1) =	- h1_5 / 120.0 
	Asys(6,2) = h2_5 / 120.0 
	Asys(6,3) =	- d6 / 720.0 
	Asys(6,4) =	h1_5 / 720.0 
	Asys(6,5) =	- h2_5 / 720.0 
	Asys(6,6) =	- n6 / 720.0 
	
	Bsys(:) = (/ -1.0, 0.0, 0.0, 0.0, 0.0, 0.0 /) 
  
    call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

	alpha = Csys(1) 
	beta  = Csys(2) 
	a = Csys(3) 
	b = Csys(4) 
	c = Csys(5) 
	d = Csys(6) 

	tri_l(k+1) = alpha 
	tri_d(k+1) = 1.0 
	tri_u(k+1) = beta 
	tri_b(k+1) = a * cell_values(offset+k-1) + b * cell_values(offset+k) + c * cell_values(offset+k+1) + d * cell_values(offset+k+2) 
 	
  end do ! end loop on cells
  
  ! Use a right-biased stencil for the second row
  
  ! Cell widths
  h0 = h(offset+1) + eps 
  h1 = h(offset+2) + eps 
  h2 = h(offset+3) + eps 
  h3 = h(offset+4) + eps 
	
  ! Auxiliary calculations
  h1_2 = h1 * h1 
  h1_3 = h1_2 * h1  
  h1_4 = h1_2 * h1_2 
  h1_5 = h1_3 * h1_2 
  h1_6 = h1_3 * h1_3 
	
  h2_2 = h2 * h2 
  h2_3 = h2_2 * h2  
  h2_4 = h2_2 * h2_2 
  h2_5 = h2_3 * h2_2 
  h2_6 = h2_3 * h2_3 
	
  g   = h0 + h1 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 

  h0ph1   = h0 + h1 
  h0ph1_2 = h0ph1 * h0ph1  
  h0ph1_3 = h0ph1_2 * h0ph1  
  h0ph1_4 = h0ph1_2 * h0ph1_2  
  h0ph1_5 = h0ph1_3 * h0ph1_2  

  d2 = ( h1_2 - g_2 ) / h0 
  d3 = ( h1_3 - g_3 ) / h0 
  d4 = ( h1_4 - g_4 ) / h0 
  d5 = ( h1_5 - g_5 ) / h0 
  d6 = ( h1_6 - g_6 ) / h0 
	
  g   = h2 + h3 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 
	
  n2 = ( g_2 - h2_2 ) / h3 
  n3 = ( g_3 - h2_3 ) / h3 
  n4 = ( g_4 - h2_4 ) / h3 
  n5 = ( g_5 - h2_5 ) / h3 
  n6 = ( g_6 - h2_6 ) / h3 
	
  ! Compute matrix entries
  Asys(1,1) = 1.0 
  Asys(1,2) = 1.0 
  Asys(1,3) = -1.0 
  Asys(1,4) = -1.0 
  Asys(1,5) = -1.0 
  Asys(1,6) = -1.0 
	
  Asys(2,1) = - h0ph1 
  Asys(2,2) = 0.0 
  Asys(2,3) = -0.5 * d2 
  Asys(2,4) = 0.5 * h1 
  Asys(2,5) = -0.5 * h2 
  Asys(2,6) = -0.5 * n2 
	
  Asys(3,1) = 0.5 * h0ph1_2 
  Asys(3,2) = 0.0 
  Asys(3,3) = d3 / 6.0 
  Asys(3,4) = - h1_2 / 6.0 
  Asys(3,5) = - h2_2 / 6.0 
  Asys(3,6) = - n3 / 6.0 
	
  Asys(4,1) = - h0ph1_3 / 6.0 
  Asys(4,2) = 0.0 
  Asys(4,3) = - d4 / 24.0 
  Asys(4,4) = h1_3 / 24.0 
  Asys(4,5) = - h2_3 / 24.0 
  Asys(4,6) = - n4 / 24.0 
	
  Asys(5,1) = h0ph1_4 / 24.0 
  Asys(5,2) = 0.0 
  Asys(5,3) = d5 / 120.0 
  Asys(5,4) = - h1_4 / 120.0 
  Asys(5,5) = - h2_4 / 120.0 
  Asys(5,6) = - n5 / 120.0 
	
  Asys(6,1) = - h0ph1_5 / 120.0 
  Asys(6,2) = 0.0 
  Asys(6,3) = - d6 / 720.0 
  Asys(6,4) = h1_5 / 720.0 
  Asys(6,5) = - h2_5 / 720.0 
  Asys(6,6) = - n6 / 720.0 
	
  Bsys(:) = (/ -1.0, h1, -0.5*h1_2, h1_3/6.0, -h1_4/24.0, h1_5/120.0 /) 

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

  alpha = Csys(1) 
  beta  = Csys(2) 
  a = Csys(3) 
  b = Csys(4) 
  c = Csys(5) 
  d = Csys(6) 

  tri_l(2) = alpha 
  tri_d(2) = 1.0 
  tri_u(2) = beta 
  tri_b(2) = a * cell_values(offset+1) + b * cell_values(offset+2) + c * cell_values(offset+3) + d * cell_values(offset+4) 
  
  ! Boundary conditions: left boundary
  x(1) = 0.0 
  do i = 2,7
  	x(i) = x(i-1) + h(offset+i-1) + eps 
  end do

  do i = 1,6
    
	do j = 1,6
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j 
	end do
	
	Bsys(i) = cell_values(offset+i) * ( h(offset+i) + eps ) 
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		
  
  tri_l(1) = 0.0 
  tri_d(1) = 1.0 
  tri_u(1) = 0.0 
  tri_b(1) = evaluation_polynomial ( Csys, 6, x(1) ) 		! first edge value
  
  ! Use a left-biased stencil for the second to last row
  
  ! Cell widths
  h0 = h(offset+N-3) + eps 
  h1 = h(offset+N-2) + eps 
  h2 = h(offset+N-1) + eps 
  h3 = h(offset+N) + eps 
	
  ! Auxiliary calculations
  h1_2 = h1 * h1 
  h1_3 = h1_2 * h1  
  h1_4 = h1_2 * h1_2 
  h1_5 = h1_3 * h1_2 
  h1_6 = h1_3 * h1_3 
	
  h2_2 = h2 * h2 
  h2_3 = h2_2 * h2  
  h2_4 = h2_2 * h2_2 
  h2_5 = h2_3 * h2_2 
  h2_6 = h2_3 * h2_3 
	
  g   = h0 + h1 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 

  h2ph3   = h2 + h3 
  h2ph3_2 = h2ph3 * h2ph3  
  h2ph3_3 = h2ph3_2 * h2ph3  
  h2ph3_4 = h2ph3_2 * h2ph3_2  
  h2ph3_5 = h2ph3_3 * h2ph3_2  

  d2 = ( h1_2 - g_2 ) / h0 
  d3 = ( h1_3 - g_3 ) / h0 
  d4 = ( h1_4 - g_4 ) / h0 
  d5 = ( h1_5 - g_5 ) / h0 
  d6 = ( h1_6 - g_6 ) / h0 
	
  g   = h2 + h3 
  g_2 = g * g 
  g_3 = g * g_2 
  g_4 = g_2 * g_2 
  g_5 = g_4 * g 
  g_6 = g_3 * g_3 
	
  n2 = ( g_2 - h2_2 ) / h3 
  n3 = ( g_3 - h2_3 ) / h3 
  n4 = ( g_4 - h2_4 ) / h3 
  n5 = ( g_5 - h2_5 ) / h3 
  n6 = ( g_6 - h2_6 ) / h3 

  ! Compute matrix entries
  Asys(1,1) = 1.0 
  Asys(1,2) = 1.0 
  Asys(1,3) = -1.0 
  Asys(1,4) = -1.0 
  Asys(1,5) = -1.0 
  Asys(1,6) = -1.0 
	
  Asys(2,1) =	0.0 
  Asys(2,2) = h2ph3 
  Asys(2,3) =	-0.5 * d2 
  Asys(2,4) =	0.5 * h1 
  Asys(2,5) =	-0.5 * h2 
  Asys(2,6) =	-0.5 * n2 
	
  Asys(3,1) =	0.0 
  Asys(3,2) = 0.5 * h2ph3_2 
  Asys(3,3) =	d3 / 6.0 
  Asys(3,4) =	- h1_2 / 6.0 
  Asys(3,5) =	- h2_2 / 6.0 
  Asys(3,6) =	- n3 / 6.0 
	
  Asys(4,1) =	0.0 
  Asys(4,2) = h2ph3_3 / 6.0 
  Asys(4,3) =	- d4 / 24.0 
  Asys(4,4) =	h1_3 / 24.0 
  Asys(4,5) =	- h2_3 / 24.0 
  Asys(4,6) =	- n4 / 24.0 
	
  Asys(5,1) =	0.0 
  Asys(5,2) = h2ph3_4 / 24.0 
  Asys(5,3) =	d5 / 120.0 
  Asys(5,4) =	- h1_4 / 120.0 
  Asys(5,5) =	- h2_4 / 120.0 
  Asys(5,6) =	- n5 / 120.0 
	
  Asys(6,1) =	0.0 
  Asys(6,2) = h2ph3_5 / 120.0 
  Asys(6,3) =	- d6 / 720.0 
  Asys(6,4) =	h1_5 / 720.0 
  Asys(6,5) =	- h2_5 / 720.0 
  Asys(6,6) =	- n6 / 720.0 

  Bsys(:) = (/ -1.0, -h2, -0.5*h2_2, -h2_3/6.0, -h2_4/24.0, -h2_5/120.0 /) 

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		

  alpha = Csys(1) 
  beta  = Csys(2) 
  a = Csys(3) 
  b = Csys(4) 
  c = Csys(5) 
  d = Csys(6) 

  tri_l(N) = alpha 
  tri_d(N) = 1.0 
  tri_u(N) = beta 
  tri_b(N) = a * cell_values(offset+N-3) + b * cell_values(offset+N-2) + c * cell_values(offset+N-1) + d * cell_values(offset+N) 

  ! Boundary conditions: right boundary
  x(1) = 0.0 
  do i = 2,7
  	x(i) = x(i-1) + h(offset+N-7+i) + eps 
  end do

  do i = 1,6
    
	do j = 1,6
	  Asys(i,j) = ( (x(i+1)**j) - (x(i)**j) ) / j 
	end do
	
	Bsys(i) = cell_values(offset+N-6+i) * ( h(offset+N-6+i) + eps ) 
	
  end do	

  call solve_linear_system ( Asys, Bsys, Csys, 6 ) 		
  
  tri_l(N+1) = 0.0 
  tri_d(N+1) = 1.0 
  tri_u(N+1) = 0.0 
  tri_b(N+1) = evaluation_polynomial ( Csys, 6, x(7) ) 		! last edge value
  
  ! Solve tridiagonal system and assign edge values
  call solve_tridiagonal_system ( tri_l, tri_d, tri_u, tri_b, tri_x, N+1 ) 

  do i = 2,N
    edge_values(offset+i,1)   = tri_x(i) 
    edge_values(offset+i-1,2) = tri_x(i) 
  end do
  edge_values(offset+1,1) = tri_x(1) 
  edge_values(offset+N,2) = tri_x(N+1)           

  deallocate(tri_l,tri_d,tri_u,tri_b,tri_x) 


end subroutine edge_values_implicit_h6

!------------------------------------------------------------------------------
! Limit pqm
! -----------------------------------------------------------------------------
subroutine pqm_limiter(nsc,nec,nsd,nsh, h, cell_values, edge_slopes,edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsd,nsh
  real,    intent(in)  :: h(nsd:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(inout) :: edge_slopes(nsh:,1:), edge_values(nsh:,1:)

!------------------------------------------------------------------------------
! Standard PQM limiter (White & Adcroft, JCP 2008).
!
! grid:  one-dimensional grid (see grid.F90)
! ppoly: piecewise quadratic polynomial to be reconstructed (see ppoly.F90)
! u: 	 cell averages
!
! It is assumed that the dimension of 'u' is equal to the number of cells
! defining 'grid' and 'ppoly'. No consistency check is performed.
!------------------------------------------------------------------------------

  ! Local variables
  integer	:: k;			! loop index
  integer	:: inflexion_l;	
  integer	:: inflexion_r;	
  real		:: u0_l, u0_r;	! edge values
  real		:: u1_l, u1_r;	! edge slopes
  real		:: u_l, u_c, u_r;		! left, center and right cell averages
  real		:: h_l, h_c, h_r;		! left, center and right cell widths
  real		:: sigma_l, sigma_c, sigma_r;	! left, center and right 
  											! van Leer slopes 	
  real		:: slope;		! retained PLM slope
  real		:: a, b, c, d, e
  real		:: alpha1, alpha2, alpha3
  real		:: rho, sqrt_rho
  real		:: gradient1, gradient2
  real		:: x1, x2
  real		:: eps;

  eps = 1e-10;
  

  ! Bound edge values
  call bound_edge_values(nsc,nec,nsd,nsh, h, cell_values, edge_values )
  ! Make discontinuous edge values monotonic (thru averaging)
  call check_discontinuous_edge_values(nsc,nec,nsh, cell_values, edge_values )

  ! Loop on interior cells to apply the PQM limiter
  do k = nsc,nec
	
	inflexion_l = 0;
	inflexion_r = 0;
	
	! Get edge values, edge slopes and cell width
	u0_l = edge_values(k,1);
	u0_r = edge_values(k,2);
	u1_l = edge_slopes(k,1);
	u1_r = edge_slopes(k,2);
	
	! Get cell widths and cell averages (boundary cells are assumed to
	! be local extrema for the sake of slopes)
	h_l = h(k-1)
	h_c = h(k) 
	h_r = h(k+1)
	u_l = cell_values(k-1);
	u_c = cell_values(k);
	u_r = cell_values(k+1);
	
	! Compute limited slope
	sigma_l = 2.0 * ( u_c - u_l ) / h_c;
	sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r );
	sigma_r = 2.0 * ( u_r - u_c ) / h_c;

	if ( (sigma_l * sigma_r) .GT. 0.0 ) then
	  slope = sign ( min (abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c );
	else
	  slope = 0.0;
	end if
	
	! If one of the slopes has the wrong sign compared with the 
	! limited PLM slope, it is set equal to the limited PLM slope
	if ( u1_l*slope .LE. 0.0 ) then
	  u1_l = slope;
	end if  
	
	if ( u1_r*slope .LE. 0.0 ) then
	  u1_r = slope;
	end if  

	! Local extremum --> flatten
	if ( (u0_r - u_c) * (u_c - u0_l) .LE. 0.0) then
	  u0_l = u_c;
	  u0_r = u_c;
	  u1_l = 0.0;
	  u1_r = 0.0;
	  inflexion_l = -1;
	  inflexion_r = -1;
	end if
	  
	! Edge values are bounded and averaged when discontinuous and not 
	! monotonic, edge slopes are consistent and the cell is not an extremum.
	! We now need to check and enforce the monotonicity of the quartic within
	! the cell
    if ( (inflexion_l .EQ. 0) .AND. (inflexion_r .EQ. 0) ) then
	
	  a = u0_l;
	  b = h_c * u1_l;
	  c = 30.0 * cell_values(k) - 12.0*u0_r - 18.0*u0_l + 1.5*h_c*(u1_r - 3.0*u1_l);
	  d = -60.0 * cell_values(k) + h_c *(6.0*u1_l - 4.0*u1_r) + 28.0*u0_r + 32.0*u0_l;
	  e = 30.0 * cell_values(k) + 2.5*h_c*(u1_r - u1_l) - 15.0*(u0_l + u0_r);
	
	  ! Determine the coefficients of the second derivative 
	  ! alpha1 xi^2 + alpha2 xi + alpha3
 	  alpha1 = 6*e
	  alpha2 = 3*d
	  alpha3 = c

	  rho = alpha2 * alpha2 - 4.0 * alpha1 * alpha3
	
	  ! Check whether inflexion points exist
	  if ( rho .GE. 0.0 ) then
	
            sqrt_rho = sqrt ( rho )
	    x1 = 0.5 * ( - alpha2 - sqrt_rho ) / (alpha1+eps)
	    x2 = 0.5 * ( - alpha2 + sqrt_rho ) / (alpha1+eps)

	    ! Check whether both inflexion points lie in [0,1]
	    if ( (x1 .GE. 0.0) .AND. (x1 .LE. 1.0) .AND. &
		     (x2 .GE. 0.0) .AND. (x2 .LE. 1.0) ) then
	  
	      gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b
	      gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b

		  ! Check whether one of the gradients is inconsistent
		  if ( (gradient1 * slope .LT. 0.0) .OR. &
		       (gradient2 * slope .LT. 0.0) ) then
		    ! Decide where to collapse inflexion points 
			! (depends on one-sided slopes)
		    if ( abs(sigma_l) .LT. abs(sigma_r) ) then
		      inflexion_l = 1
		    else
		      inflexion_r = 1
		    end if
		  end if
	
	    ! If both x1 and x2 do not lie in [0,1], check whether 
		! only x1 lies in [0,1]
	    else if ( (x1 .GE. 0.0) .AND. (x1 .LE. 1.0) ) then
	 
	      gradient1 = 4.0 * e * (x1**3) + 3.0 * d * (x1**2) + 2.0 * c * x1 + b

	 	  ! Check whether the gradient is inconsistent
		  if ( gradient1 * slope .LT. 0.0 ) then
		    ! Decide where to collapse inflexion points 
			! (depends on one-sided slopes)
		    if ( abs(sigma_l) .LT. abs(sigma_r) ) then
		      inflexion_l = 1
		    else
		      inflexion_r = 1
		    end if
		  end if
	
	    ! If x1 does not lie in [0,1], check whether x2 lies in [0,1]
	    else if ( (x2 .GE. 0.0) .AND. (x2 .LE. 1.0) ) then
	    
		  gradient2 = 4.0 * e * (x2**3) + 3.0 * d * (x2**2) + 2.0 * c * x2 + b
		
		  ! Check whether the gradient is inconsistent
		  if ( gradient2 * slope .LT. 0.0 ) then
		    ! Decide where to collapse inflexion points 
			! (depends on one-sided slopes)
		    if ( abs(sigma_l) .LT. abs(sigma_r) ) then
		      inflexion_l = 1
		    else
		      inflexion_r = 1
		    end if
		  end if

	    end if ! end checking where the inflexion points lie 

	  end if ! end checking if rho >= 0
	   
	end if ! end checking whether to shift inflexion points

	! At this point, we know onto which edge to shift inflexion points
	if ( inflexion_l .EQ. 1 ) then
	  
	  ! We modify the edge slopes so that both inflexion points 
	  ! collapse onto the left edge
	  u1_l = ( 10.0 * u_c - 2.0 * u0_r - 8.0 * u0_l ) / (3.0*h_c)
	  u1_r = ( -10.0 * u_c + 6.0 * u0_r + 4.0 * u0_l ) / h_c

	  ! One of the modified slopes might be inconsistent. When that happens,
	  ! the inconsistent slope is set equal to zero and the opposite edge value
	  ! and edge slope are modified in compliance with the fact that both 
	  ! inflexion points must still be located on the left edge
	  if ( u1_l * slope .LT. 0.0 ) then
	  
	    u1_l = 0.0
		u0_r = 5.0 * u_c - 4.0 * u0_l
		u1_r = 20.0 * (u_c - u0_l) / h_c
	  
	  else if ( u1_r * slope .LT. 0.0 ) then
	 	
		u1_r = 0.0
		u0_l = (5.0*u_c - 3.0*u0_r) / 2.0
		u1_l = 10.0 * (-u_c + u0_r) / (3.0 * h_c)
	  
	  end if
	  
	else if ( inflexion_r .EQ. 1 ) then
	  
	  ! We modify the edge slopes so that both inflexion points 
	  ! collapse onto the right edge
	  u1_r = ( -10.0 * u_c + 8.0 * u0_r + 2.0 * u0_l ) / (3.0 * h_c)
	  u1_l = ( 10.0 * u_c - 4.0 * u0_r - 6.0 * u0_l ) / h_c
	  
	  ! One of the modified slopes might be inconsistent. When that happens,
	  ! the inconsistent slope is set equal to zero and the opposite edge value
	  ! and edge slope are modified in compliance with the fact that both 
	  ! inflexion points must still be located on the right edge
	  if ( u1_l * slope .LT. 0.0 ) then
	    
		u1_l = 0.0
		u0_r = ( 5.0 * u_c - 3.0 * u0_l ) / 2.0
		u1_r = 10.0 * (u_c - u0_l) / (3.0 * h_c)
	  
	  else if ( u1_r * slope .LT. 0.0 ) then
	    
		u1_r = 0.0
		u0_l = 5.0 * u_c - 4.0 * u0_r
		u1_l = 20.0 * ( -u_c + u0_r ) / h_c
	  
	  end if

	end if ! clause to check where to collapse inflexion points
	
    ! Save edge values and edge slopes for reconstruction
	edge_values(k,1) = u0_l
	edge_values(k,2) = u0_r
	edge_slopes(k,1) = u1_l
	edge_slopes(k,2) = u1_r
    
  end do ! end loop on interior cells

  ! Constant reconstruction within boundary cells
  edge_values(nsc,:) = cell_values(nsc);
  edge_slopes(nsc,:) = 0.0;
  
  edge_values(nec,:) = cell_values(nec);
  edge_slopes(nec,:) = 0.0;
  
end subroutine pqm_limiter

!------------------------------------------------------------------------------
! Bound edge values by neighboring cell averages
!------------------------------------------------------------------------------
subroutine bound_edge_values (nsc,nec,nsd,nsh, h, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsd,nsh
  real,    intent(in)  :: h(nsd:) 	! cell width
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)

! ------------------------------------------------------------------------------
! In this routine, we loop on all cells to bound their left and right
! edge values by the cell averages. That is, the left edge value must lie
! between the left cell average and the central cell average. A similar
! reasoning applies to the right edge values. 
!
! Both boundary edge values are set equal to the boundary cell averages.
! Any extrapolation scheme is applied after this routine has been called.
! Therefore, boundary cells are treated as if they were local extrema.
! ------------------------------------------------------------------------------

  ! Local variables
  integer		:: k;			! loop index
  real			:: h_l, h_c, h_r;
  real			:: u_l, u_c, u_r;
  real			:: u0_l, u0_r;
  real			:: sigma_l, sigma_c, sigma_r;	! left, center and right 
  												! van Leer slopes 	
  real			:: slope;				! retained PLM slope

  ! Loop on cells to bound edge value
!  do k = nsc+1,nec-1
  do k = nsc,nec
    
	! For the sake of bounding boundary edge values, the left neighbor
	! of the left boundary cell is assumed to be the same as the left
	! boundary cell and the right neighbor of the right boundary cell
	! is assumed to be the same as the right boundary cell. This 
	! effectively makes boundary cells look like extrema.

    ! All cells can now be treated equally
	h_l = h(k-1);
	h_c = h(k);
	h_r = h(k+1);
	
	u_l = cell_values(k-1);
	u_c = cell_values(k);
	u_r = cell_values(k+1);

	u0_l = edge_values(k,1);
	u0_r = edge_values(k,2);

	sigma_l = 2.0 * ( u_c - u_l ) / h_c;
	sigma_c = 2.0 * ( u_r - u_l ) / ( h_l + 2.0*h_c + h_r );
	sigma_r = 2.0 * ( u_r - u_c ) / h_c;
	
	if ( (sigma_l * sigma_r) .GT. 0.0 ) then
	  slope = sign ( min (abs(sigma_l),abs(sigma_c),abs(sigma_r)), sigma_c );
	else
	  slope = 0.0;
	end if
	
	! The limiter must be used in the local coordinate system to each cell.
	! Hence, we must multiply the slope by h1. The multiplication by 0.5 is
	! simply a way to make it useable in the limiter (cfr White and Adcroft 
	! JCP 2008 Eqs 19 and 20)
	slope = slope * h_c * 0.5;
	
	if ( (u_l-u0_l)*(u0_l-u_c) .LT. 0.0 ) then
	  u0_l = u_c - sign ( min ( abs(slope), abs(u0_l-u_c) ), slope );
	end if
	
	if ( (u_r-u0_r)*(u0_r-u_c) .LT. 0.0 ) then
	  u0_r = u_c + sign ( min ( abs(slope), abs(u0_r-u_c) ), slope );
	end if

	! Store edge values
	edge_values(k,1) = u0_l;
	edge_values(k,2) = u0_r;

  end do ! loop on interior edges

end subroutine bound_edge_values


!------------------------------------------------------------------------------
! Check discontinuous edge values and take average is not monotonic
!------------------------------------------------------------------------------
subroutine check_discontinuous_edge_values (nsc,nec,nsh, cell_values, edge_values )
  ! Arguments
  integer, intent(in)  :: nsc,nec,nsh
  real,    intent(in)  :: cell_values(nsh:)
  real,    intent(out) :: edge_values(nsh:,1:)
! ------------------------------------------------------------------------------
! For each interior edge, check whether the edge values are discontinuous.
! If so and if they are not monotonic, replace each edge value by their average.
! ------------------------------------------------------------------------------

  ! Local variables
  integer		:: k;			! loop index
  real			:: u0_minus;	! left value at given edge
  real			:: u0_plus;		! right value at given edge
  real			:: um_minus;	! left cell average
  real			:: um_plus;		! right cell average
  real			:: u0_avg;		! avg value at given edge


  ! Loop on interior cells
!  do k = nsc,nec-1
  do k = nsc-1,nec
  
    ! Edge value on the left of the edge
	u0_minus = edge_values(k,2);

	! Edge value on the right of the edge
	u0_plus  = edge_values(k+1,1);

	! Left cell average
	um_minus = cell_values(k);

	! Right cell average
	um_plus  = cell_values(k+1);

	if ( (u0_plus - u0_minus)*(um_plus - um_minus) .LT. 0.0 ) then
	  u0_avg = 0.5 * ( u0_minus + u0_plus );
	  edge_values(k  ,2) = u0_avg;
	  edge_values(k+1,1) = u0_avg;
	end if
	
  end do ! end loop on interior edges

end subroutine check_discontinuous_edge_values

! -----------------------------------------------------------------------------
! Pointwise evaluation of a polynomial
! -----------------------------------------------------------------------------
real function evaluation_polynomial ( coefficients, nb_coefficients, x )
! -----------------------------------------------------------------------------
! The polynomial is defined by the coefficients contained in the 
! array of the same name, as follows: C(1) + C(2)x + C(3)x^2 + C(4)x^3 + ... 
! where C refers to the array 'coefficients'.
! The number of coefficients is given by nb_coefficients and x
! is the coordinate where the polynomial is to be evaluated.
!
! The function returns the value of the polynomial at x.
! -----------------------------------------------------------------------------
  
  ! Arguments
  real, dimension(:), intent(in)		:: coefficients
  integer, intent(in)					:: nb_coefficients
  real, intent(in)						:: x
  
  ! Local variables
  integer								:: k	
  real									:: f 	! value of polynomial at x
  
  f = 0.0
  do k = 1,nb_coefficients
  	f = f + coefficients(k) * ( x**(k-1) )
  end do

  evaluation_polynomial = f

end function evaluation_polynomial


! -----------------------------------------------------------------------------
! Exact integration of polynomial of degree n
! -----------------------------------------------------------------------------
real function integration_polynomial ( xi0, xi1, C, n )
! -----------------------------------------------------------------------------
! Exact integration of a polynomial of degree n over the interval [xi0,xi1].
! The array of coefficients (C) must be of size n+1, where n is the degree of
! the polynomial to integrate.
! -----------------------------------------------------------------------------
  
  ! Arguments
  real, intent(in)					:: xi0, xi1
  real, dimension(:), intent(in)	:: C
  integer, intent(in)				:: n

  ! Local variables
  integer							:: k	
  real								:: integral
  
  integral = 0.0
  
  do k = 1,(n+1) 
    integral = integral + C(k) * (xi1**k - xi0**k) / real(k);
  end do

  integration_polynomial = integral
  
end function integration_polynomial

! -----------------------------------------------------------------------------
! Solve the linear system AX = B
! -----------------------------------------------------------------------------
subroutine solve_linear_system ( A, B, X, system_size )
! -----------------------------------------------------------------------------
! This routine uses Gauss's algorithm to transform the system's original 
! matrix into an upper triangular matrix. Back substitution yields the answer.
! The matrix A must be square and its size must be that of the vectors B and X.
! -----------------------------------------------------------------------------

  ! Arguments
  real, dimension(:,:), intent(inout)	:: A;
  real, dimension(:), intent(inout)		:: B;
  real, dimension(:), intent(inout)		:: X;
  integer								:: system_size;

  ! Local variables
  integer				:: i, j, k;
  real, parameter		:: eps = 0.0; 		! Minimum pivot magnitude allowed
  real					:: factor;
  real					:: pivot;
  real					:: swap_a, swap_b;
  integer				:: found_pivot;		! boolean indicating whether
  											! a pivot has been found
  ! Loop on rows
  do i = 1,system_size-1

    found_pivot = 0;

    ! Start to look for a pivot in row i. If the pivot
	! in row i -- which is the current row -- is not valid,
	! we keep looking for a valid pivot by searching the
	! entries of column i in rows below row i. Once a valid 
	! pivot is found (say in row k), rows i and k are swaped.
	k = i;	
	do while ( ( found_pivot .ne. 0 ) .AND. ( k .LE. system_size ) )
	    
		if ( abs( A(k,i) ) .GT. eps ) then	! a valid pivot is found
		  found_pivot = 1;
		else								! Go to the next row to see 
											! if there is a valid pivot there
		  k = k + 1;	
		end if
		
	end do ! end loop to find pivot

	! If no pivot could be found, the system is singular and we need
	! to end the execution
	if ( found_pivot .ne. 0 ) then
	  write(*,*) 'The linear system is singular !'
          exit
	end if

	! If the pivot is in a row that is different than row i, that is if
	! k is different than i, we need to swap those two rows
	if ( k .NE. i ) then
	  do j = 1,system_size
	    swap_a = A(i,j);
	    A(i,j) = A(k,j);
	    A(k,j) = swap_a;
	  end do
	  swap_b = B(i);
	  B(i) = B(k);
	  B(k) = swap_b;	
	end if
			
	! Trasnform pivot to 1 by dividing the entire row 
	! (right-hand side included) by the pivot
	pivot = A(i,i)

        if(abs(pivot).lt.1.0e-10) write(6,*) 'pivot= ',pivot

	do j = i,system_size
	  A(i,j) = A(i,j) / pivot;
	end do
	B(i) = B(i) / pivot;

	! #INV: At this point, A(i,i) is a suitable pivot and it is equal to 1
			
	! Put zeros in column for all rows below that containing 
	! pivot (which is row i)
	do k = (i+1),system_size	! k is the row index
	  factor = A(k,i);
	  do j = (i+1),system_size		! j is the column index
	    A(k,j) = A(k,j) - factor * A(i,j);
	  end do
	  B(k) = B(k) - factor * B(i);
	end do
  
  end do ! end loop on i	
	

  ! Solve system by back substituting
  X(system_size) = B(system_size) / A(system_size,system_size);
  do i = system_size-1,1,-1	! loop on rows, starting from second to last row
    X(i) = B(i);
	do j = (i+1),system_size 
	  X(i) = X(i) - A(i,j) * X(j);
	end do
	X(i) = X(i) / A(i,i);
  end do

end subroutine solve_linear_system


! -----------------------------------------------------------------------------
! Solve the tridiagonal system AX = B 
! -----------------------------------------------------------------------------
subroutine solve_tridiagonal_system ( Al, Ad, Au, B, X, system_size )
! -----------------------------------------------------------------------------
! This routine uses Thomas's algorithm to solve the tridiagonal system AX = B.
! (A is made up of lower, middle and upper diagonals)
! -----------------------------------------------------------------------------
  ! Arguments	
  real, dimension(:), intent(inout)	:: Al, Ad, Au;	! lo., mid. and up. diagonals
  real, dimension(:), intent(inout)	:: B;			! system right-hand side
  real, dimension(:), intent(inout)	:: X;			! solution vector	
  integer, intent(in)				:: system_size;
	
  ! Local variables	
  integer								:: k;	    ! Loop index
  integer								:: N;		! system size 

  N = system_size;	

  ! Factorization
  do k = 1,N-1
	Al(k+1) = Al(k+1) / Ad(k);
	Ad(k+1) = Ad(k+1) - Al(k+1) * Au(k);
  end do

  ! Forward sweep
  do k = 2,N
  	B(k) = B(k) - Al(k) * B(k-1);
  end do

  ! Backward sweep
  X(N) = B(N) / Ad(N);
  do k = N-1,1,-1
    X(k) = ( B(k) - Au(k)*X(k+1) ) / Ad(k);
  end do

end subroutine solve_tridiagonal_system

subroutine xTracerFlx_OS7MP_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_OS7MP_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real    :: Eps,Fac,DelP,Phi,DelM,Del2,DelPP,Del2P,Del3P,Del4,DelPPP,Del2PP,Del3PP,Del2M,Del3M,DelMMM
  real    :: Del4P,Del5P,DelMM,Del2MM,Del3MM,Del4M,Del5M,Del6,DelIP,DelI,recip_DelIp,recip_DelI
  real    :: rp1h,rp1h_cfl,d2,d2p1,d2m1,A,B,C,D,dp1h,dm1h,PhiMD,PhiLC,PhiMin,PhiMax,Psi
  real    :: Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm,MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
  real	  :: cfl, massflux

  tracer_flux = 0.0

  Eps = 1.0d-20

  do j=JNLO_,JNHI_ ; do I=IS_,IE_ 

     massflux = dy(I,j) * uhrho(I,j)

     cfl = abs( u(I,j) * dt /dx(I,j)) 

     if (massflux.gt.0.0) then                              
        Qippp = tracer_field(i+3,j)
        Qipp  = tracer_field(i+2,j)
        Qip   = tracer_field(i+1,j)
        Qi    = tracer_field(i,j)
        Qim   = tracer_field(i-1,j)
        Qimm  = tracer_field(i-2,j)
        Qimmm = tracer_field(i-3,j)

        MskIpp  = tmask(i+2,j)
        MskIp   = tmask(i+1,j)
        MskI    = tmask(i,j)
        MskIm   = tmask(i-1,j)
        MskImm  = tmask(i-2,j)
        MskImmm = tmask(i-3,j)
     elseif (massflux.lt.0.0) then
        Qippp = tracer_field(i-2,j)
        Qipp  = tracer_field(i-1,j)
        Qip   = tracer_field(i,j)
        Qi    = tracer_field(i+1,j)
        Qim   = tracer_field(i+2,j)
        Qimm  = tracer_field(i+3,j)
        Qimmm = tracer_field(i+4,j)

        MskIpp  = tmask(i-1,j)
        MskIp   = tmask(i,j)
        MskI    = tmask(i+1,j)
        MskIm   = tmask(i+2,j)
        MskImm  = tmask(i+3,j)
        MskImmm = tmask(i+4,j)
     else
        Qippp = 0.
        Qipp  = 0.
        Qip   = 0.
        Qi    = 0.
        Qim   = 0.
        Qimm  = 0.
        Qimmm = 0.

        MskIpp  = 0.
        MskIp   = 0.
        MskI    = 0.
        MskIm   = 0.
        MskImm  = 0.
        MskImmm = 0.
     endif

     !        2nd order correction [i i-1]
     Fac = 1. 
     DelP = (Qip-Qi)*MskI
     Phi = Fac * DelP
     !        3rd order correction [i i-1 i-2]
     Fac = Fac * ( cfl + 1.  )/3. 
     DelM = (Qi-Qim)*MskIm
     Del2 = DelP - DelM
     Phi = Phi - Fac * Del2
     !        4th order correction [i+1 i i-1 i-2]
     Fac = Fac * ( cfl - 2.  )/4. 
     DelPP = (Qipp-Qip)*MskIp*MskI
     Del2P = DelPP - DelP
     Del3P = Del2P - Del2
     Phi = Phi + Fac * Del3p
     !        5th order correction [i+1 i i-1 i-2 i-3]
     Fac = Fac * ( cfl - 3.  )/5. 
     DelMM = (Qim-Qimm)*MskImm*MskIm
     Del2M = DelM - DelMM
     Del3M = Del2 - Del2M
     Del4 = Del3P - Del3M
     Phi = Phi + Fac * Del4
     !        6th order correction [i+2 i+1 i i-1 i-2 i-3]
     Fac = Fac * ( cfl + 2.  )/6. 
     DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
     Del2PP = DelPP - DelP
     Del3PP = Del2PP - Del2P
     Del4P = Del3PP - Del3P
     Del5P = Del4P - Del4
     Phi = Phi + Fac * Del5P
     !        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
     Fac = Fac * ( cfl + 2.  )/7. 
     DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
     Del2MM = DelMM - DelMMM
     Del3MM = Del2M - Del2MM
     Del4M = Del3M - Del3MM
     Del5M = Del4 - Del4M
     Del6 = Del5P - Del5M
     Phi = Phi - Fac * Del6

     DelIp = ( Qip - Qi ) * MskI
     !        Phi = sign(1. ,Phi)*sign(1. ,DelIp)
     !    &        *abs(Phi+Eps)/abs(DelIp+Eps)
     !--   simplify and avoid division by zero
     recip_DelIp = sign(1. ,DelIp)/max(abs(DelIp),Eps)
     Phi = Phi*recip_DelIp

     DelI = ( Qi - Qim ) * MskIm
     !        rp1h =sign(1. ,DelI)*sign(1. ,DelIp)
     !    &        *abs(DelI+Eps)/abs(DelIp+Eps)
     !--   simplify and avoid division by zero
     recip_DelI = sign(1. ,DelI)/max(abs(DelI),Eps)
     rp1h = DelI*recip_DelIp
     rp1h_cfl = rp1h/(cfl+Eps)

     !        TVD limiter
     !        Phi = max(0. , min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )

     !        MP limiter
     d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
     d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
     d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
     A = 4. *d2 - d2p1
     B = 4. *d2p1 - d2
     C = d2
     D = d2p1
     dp1h = max(min(A,B,C,D),0. )+min(max(A,B,C,D),0. )
     A = 4. *d2m1 - d2
     B = 4. *d2 - d2m1
     C = d2m1
     D = d2
     dm1h = max(min(A,B,C,D),0. )+min(max(A,B,C,D),0. )
     !        qMD = 0.5*( ( Qi + Qip ) - dp1h )
     !        qMD = 0.5 *( ( 2. *Qi + DelIp ) - dp1h )
     !        qUL = Qi + (1. -cfl)/(cfl+Eps)*DelI
     !        qLC = Qi + 0.5 *( 1. +dm1h/(DelI+Eps) )*(qUL-Qi)
     !        PhiMD = 2. /(1. -cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
     !        PhiLC = 2. *rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
     !--   simplify and avoid division by zero
     PhiMD = 1. /(1. -cfl)*(DelIp-dp1h)*recip_DelIp
     PhiLC = rp1h_cfl*( 1. +dm1h*recip_DelI )
     !--
     PhiMin = max( min(0. ,PhiMD), min(0. ,2. *rp1h_cfl,PhiLC) )
     PhiMax = min( max(2. /(1. -cfl),PhiMD), max(0. ,2. *rp1h_cfl,PhiLC) )
     Phi = max(PhiMin,min(Phi,PhiMax))

     Psi = Phi * 0.5  * (1.  - cfl)
     tracer_flux(I,j) = massflux*( Qi + Psi*DelIp )

  enddo;enddo

end subroutine xTracerFlx_OS7MP_cne

subroutine yTracerFlx_OS7MP_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: yTracerFlx_OS7MP_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 4 points on each boundary.

  !Reference: This kernel was ported from MITGCM.

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of cell at u-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of face at t-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real    :: Eps,Fac,DelP,Phi,DelM,Del2,DelPP,Del2P,Del3P,Del4,DelPPP,Del2PP,Del3PP,Del2M,Del3M,DelMMM
  real    :: Del4P,Del5P,DelMM,Del2MM,Del3MM,Del4M,Del5M,Del6,DelIP,DelI,recip_DelIp,recip_DelI
  real    :: rp1h,rp1h_cfl,d2,d2p1,d2m1,A,B,C,D,dp1h,dm1h,PhiMD,PhiLC,PhiMin,PhiMax,Psi
  real    :: Qippp,Qipp,Qip,Qi,Qim,Qimm,Qimmm,MskIpp,MskIp,MskI,MskIm,MskImm,MskImmm
  real	  :: cfl, massflux

  tracer_flux = 0.0

  Eps = 1.0d-20

  do J=JS_,JE_; do i=INLO_,INHI_

     massflux = dx(i,J) * vhrho(i,J)

     cfl = abs( v(i,J) * dt /dy(i,J)) 

     if (massflux.gt.0.0) then                              
        Qippp = tracer_field(i,j+3)
        Qipp  = tracer_field(i,j+2)
        Qip   = tracer_field(i,j+1)
        Qi    = tracer_field(i,j)
        Qim   = tracer_field(i,j-1)
        Qimm  = tracer_field(i,j-2)
        Qimmm = tracer_field(i,j-3)

        MskIpp  = tmask(i,j+2)
        MskIp   = tmask(i,j+1)
        MskI    = tmask(i,j)
        MskIm   = tmask(i,j-1)
        MskImm  = tmask(i,j-2)
        MskImmm = tmask(i,j-3)
        !        ELSEIF (uTrans(i,j).LT.0. _d 0) THEN
     elseif (massflux.lt.0.0) then
        Qippp = tracer_field(i,j-2)
        Qipp  = tracer_field(i,j-1)
        Qip   = tracer_field(i,j)
        Qi    = tracer_field(i,j+1)
        Qim   = tracer_field(i,j+2)
        Qimm  = tracer_field(i,j+3)
        Qimmm = tracer_field(i,j+4)

        MskIpp  = tmask(i,j-1)
        MskIp   = tmask(i,j)
        MskI    = tmask(i,j+1)
        MskIm   = tmask(i,j+2)
        MskImm  = tmask(i,j+3)
        MskImmm = tmask(i,j+4)
     else
        Qippp = 0.
        Qipp  = 0.
        Qip   = 0.
        Qi    = 0.
        Qim   = 0.
        Qimm  = 0.
        Qimmm = 0.

        MskIpp  = 0.
        MskIp   = 0.
        MskI    = 0.
        MskIm   = 0.
        MskImm  = 0.
        MskImmm = 0.
     endif

     !        2nd order correction [i i-1]
     Fac = 1. 
     DelP = (Qip-Qi)*MskI
     Phi = Fac * DelP
     !        3rd order correction [i i-1 i-2]
     Fac = Fac * ( cfl + 1.  )/3. 
     DelM = (Qi-Qim)*MskIm
     Del2 = DelP - DelM
     Phi = Phi - Fac * Del2
     !        4th order correction [i+1 i i-1 i-2]
     Fac = Fac * ( cfl - 2.  )/4. 
     DelPP = (Qipp-Qip)*MskIp*MskI
     Del2P = DelPP - DelP
     Del3P = Del2P - Del2
     Phi = Phi + Fac * Del3p
     !        5th order correction [i+1 i i-1 i-2 i-3]
     Fac = Fac * ( cfl - 3.  )/5. 
     DelMM = (Qim-Qimm)*MskImm*MskIm
     Del2M = DelM - DelMM
     Del3M = Del2 - Del2M
     Del4 = Del3P - Del3M
     Phi = Phi + Fac * Del4
     !        6th order correction [i+2 i+1 i i-1 i-2 i-3]
     Fac = Fac * ( cfl + 2.  )/6. 
     DelPPP = (Qippp-Qipp)*MskIpp*MskIp*MskI
     Del2PP = DelPP - DelP
     Del3PP = Del2PP - Del2P
     Del4P = Del3PP - Del3P
     Del5P = Del4P - Del4
     Phi = Phi + Fac * Del5P
     !        7th order correction [i+2 i+1 i i-1 i-2 i-3 i-4]
     Fac = Fac * ( cfl + 2.  )/7. 
     DelMMM = (Qimm-Qimmm)*MskImmm*MskImm*MskIm
     Del2MM = DelMM - DelMMM
     Del3MM = Del2M - Del2MM
     Del4M = Del3M - Del3MM
     Del5M = Del4 - Del4M
     Del6 = Del5P - Del5M
     Phi = Phi - Fac * Del6

     DelIp = ( Qip - Qi ) * MskI
     !        Phi = sign(1. ,Phi)*sign(1. ,DelIp)
     !    &        *abs(Phi+Eps)/abs(DelIp+Eps)
     !--   simplify and avoid division by zero
     recip_DelIp = sign(1. ,DelIp)/max(abs(DelIp),Eps)
     Phi = Phi*recip_DelIp

     DelI = ( Qi - Qim ) * MskIm
     !        rp1h =sign(1. ,DelI)*sign(1. ,DelIp)
     !    &        *abs(DelI+Eps)/abs(DelIp+Eps)
     !--   simplify and avoid division by zero
     recip_DelI = sign(1. ,DelI)/max(abs(DelI),Eps)
     rp1h = DelI*recip_DelIp
     rp1h_cfl = rp1h/(cfl+Eps)

     !        TVD limiter
     !        Phi = max(0. , min( 2./(1-cfl), Phi, 2.*rp1h_cfl ) )

     !        MP limiter
     d2   = Del2 !( ( Qip + Qim ) - 2.*Qi  ) * MskI * MskIm
     d2p1 = Del2P !( ( Qipp + Qi ) - 2.*Qip ) * MskIp * MskI
     d2m1 = Del2M !( ( Qi + Qimm ) - 2.*Qim ) * MskIm * MskImm
     A = 4. *d2 - d2p1
     B = 4. *d2p1 - d2
     C = d2
     D = d2p1
     dp1h = max(min(A,B,C,D),0. )+min(max(A,B,C,D),0. )
     A = 4. *d2m1 - d2
     B = 4. *d2 - d2m1
     C = d2m1
     D = d2
     dm1h = max(min(A,B,C,D),0. )+min(max(A,B,C,D),0. )
     !        qMD = 0.5*( ( Qi + Qip ) - dp1h )
     !        qMD = 0.5 *( ( 2. *Qi + DelIp ) - dp1h )
     !        qUL = Qi + (1. -cfl)/(cfl+Eps)*DelI
     !        qLC = Qi + 0.5 *( 1. +dm1h/(DelI+Eps) )*(qUL-Qi)
     !        PhiMD = 2. /(1. -cfl)*(qMD-Qi+Eps)/(DelIp+Eps)
     !        PhiLC = 2. *rp1h_cfl*(qLC-Qi+Eps)/(qUL-Qi+Eps)
     !--   simplify and avoid division by zero
     PhiMD = 1. /(1. -cfl)*(DelIp-dp1h)*recip_DelIp
     PhiLC = rp1h_cfl*( 1. +dm1h*recip_DelI )
     !--
     PhiMin = max( min(0. ,PhiMD), min(0. ,2. *rp1h_cfl,PhiLC) )
     PhiMax = min( max(2. /(1. -cfl),PhiMD), max(0. ,2. *rp1h_cfl,PhiLC) )
     Phi = max(PhiMin,min(Phi,PhiMax))

     Psi = Phi * 0.5  * (1.  - cfl)
     tracer_flux(i,J) = massflux*( Qi + Psi*DelIp )

  enddo;enddo

end subroutine yTracerFlx_OS7MP_cne

subroutine xTracerFlx_DST3FL_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_DST3FL_cne(sv,tmask,dt,dx,dy,u,uhrho,tracer_field, tracer_flux)

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,JN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(UN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     u(UN_,JN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: uhrho(UN_,JN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(UN_,JN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real :: massflux, cfl
  real :: Rjm,Rj,Rjp,d0,d1,psiP,psiM,thetaP,thetaM,Eps
  real, parameter :: thetaMax=1.e+20

  Eps = 1.0e-20

  tracer_flux = 0.0

  do j=JNLO_,JNHI_ ; do I=IS_,IE_

     massflux = dy(I,j) * uhrho(I,j)

     cfl = abs( u(I,j) * dt /dx(I,j)) 

     Rjp=(tracer_field(i+2,j)-tracer_field(i+1,j))*tmask(i+2,j) *tmask(i+1,j)
     Rj =(tracer_field(i+1,j)-tracer_field(i,j)  )*tmask( i+1,j)*tmask( i ,j)
     Rjm=(tracer_field(i,j)  -tracer_field(i-1,j))*tmask(i,j)   *tmask(i-1,j)

     d0=(2.-cfl)*(1.-cfl)/6
     d1=(1.-cfl*cfl)/6

     !      the old version: can produce overflow, division by zero,
     !       and is wrong for tracer with low concentration:
     !       thetaP=Rjm/(1.D-20+Rj)
     !       thetaM=Rjp/(1.D-20+Rj)
     !      the right expression, but not bounded:
     !       thetaP=0.D0
     !       thetaM=0.D0
     !       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
     !       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
     !      prevent |thetaP,M| to reach too big value:
  
     IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
        thetaP=SIGN(thetaMax,Rjm*Rj)
     ELSE
        thetaP=Rjm/Rj
     ENDIF
     IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
        thetaM=SIGN(thetaMax,Rjp*Rj)
     ELSE
        thetaM=Rjp/Rj
     ENDIF

     psiP=d0+d1*thetaP
     psiP=MAX(0.,MIN(MIN(1.,psiP),thetaP*(1.-cfl)/(cfl+Eps) ))
     psiM=d0+d1*thetaM
     psiM=MAX(0.,MIN(MIN(1.,psiM),thetaM*(1.-cfl)/(cfl+Eps) ))

     tracer_flux(I,j)=                         &
          0.5*(massflux+ABS(massflux))           &
          *( tracer_field(i,j)  + psiP*Rj)  &
          +0.5*(massflux-ABS(massflux))           &
          *( tracer_field(i+1,j)- psiM*Rj)

  enddo;enddo

end subroutine xTracerFlx_DST3FL_cne

subroutine yTracerFlx_DST3FL_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !INTERFACE: xTracerFlx_DST3FL_cne(sv,tmask,dt,dx,dy,v,vhrho,tracer_field, tracer_flux)

  !NOTES: Input arrays tmask and tracer_field must be properly haloed by at least 2 points on each boundary.

  !Reference:

  !ARGUMENTS:
  integer, intent(in) :: SV_            	! Shape vector
  real,    intent(in) ::    dt                  ! time step increment
  real,    intent(in) ::    dx(IN_,VN_) 	! x width of cell at t-point (m)
  real,    intent(in) ::    dy(IN_,JN_) 	! y width of face at u-point (m)
  real,    intent(in) ::     v(IN_,VN_) 	! Component of velocity normal to face (m/s)
  real,    intent(in) :: vhrho(IN_,VN_) 	! Component of mass flow normal to face (Kg/m/s)
  real,    intent(in) ::        tmask(INH_,JNH_)! tracer mask
  real,    intent(in) :: tracer_field(INH_,JNH_)! Tracer field at c-points ([tracer_concentration]) 
  real,    intent(out)::  tracer_flux(IN_,VN_)! Tracer flux  at u-points ([tracer_concentration]*Kg/s)
  !LOCAL VARIABLES:
  integer :: i,j
  real :: massflux, cfl
  real :: Rjm,Rj,Rjp,d0,d1,psiP,psiM,thetaP,thetaM,Eps
  real, parameter :: thetaMax=1.d+20

  Eps = 1.0d-20
  tracer_flux = 0.0

  do J=JS_,JE_ ; do i=INLO_,INHI_

     massflux = dx(i,J) * vhrho(i,J)

     cfl = abs( v(i,J) * dt /dy(i,J)) 

     Rjp=(tracer_field(i,j+2)-tracer_field(i,j+1))*tmask(i,j+2)*tmask(i,j+1)
     Rj =(tracer_field(i,j+1)-tracer_field(i,j)  )*tmask(i,j+1)*tmask( i ,j)
     Rjm=(tracer_field(i,j)  -tracer_field(i,j-1))*tmask(i,j)  *tmask(i,j-1)

     d0=(2.-cfl)*(1.-cfl)/6
     d1=(1.-cfl*cfl)/6

     !      the old version: can produce overflow, division by zero,
     !       and is wrong for tracer with low concentration:
     !       thetaP=Rjm/(1.D-20+Rj)
     !       thetaM=Rjp/(1.D-20+Rj)
     !      the right expression, but not bounded:
     !       thetaP=0.D0
     !       thetaM=0.D0
     !       IF (Rj.NE.0.D0) thetaP=Rjm/Rj
     !       IF (Rj.NE.0.D0) thetaM=Rjp/Rj
     !      prevent |thetaP,M| to reach too big value:
  
     IF ( ABS(Rj)*thetaMax .LE. ABS(Rjm) ) THEN
        thetaP=SIGN(thetaMax,Rjm*Rj)
     ELSE
        thetaP=Rjm/Rj
     ENDIF
     IF ( ABS(Rj)*thetaMax .LE. ABS(Rjp) ) THEN
        thetaM=SIGN(thetaMax,Rjp*Rj)
     ELSE
        thetaM=Rjp/Rj
     ENDIF

     psiP=d0+d1*thetaP
     psiP=MAX(0.,MIN(MIN(1.,psiP),thetaP*(1.-cfl)/(cfl+Eps) ))
     psiM=d0+d1*thetaM
     psiM=MAX(0.,MIN(MIN(1.,psiM),thetaM*(1.-cfl)/(cfl+Eps) ))

     tracer_flux(i,J)=                         &
          0.5*(massflux+ABS(massflux))           &
          *( tracer_field(i,j)  + psiP*Rj)  &
          +0.5*(massflux-ABS(massflux))           &
          *( tracer_field(i,j+1)- psiM*Rj)

  enddo;enddo


end subroutine yTracerFlx_DST3FL_cne


end module HOME_kernels_cne
!EOP

